share/colnew/prob4.mac

   1 /* Problem 4 for colnew, example of continuation,
   2    written by Michel Talon
   3
   4    If everything is working correctly with Maxima's interface to
   5    colnew, then this should match the Fortran result in the ex4
   6    directory, which is a translation of this to Fortran.
   7 */
   8
   9 /******************************************************************/
  10 /* Quantum Neumann equation, with 2 "eigenvalues" c_1 and c_2     */
  11 /* (c_1=v-c_2-c_3, v is a parameter, used in continuation)        */
  12 /*                                                                */
  13 /* diff(f,x,2) + (1/2)*(1/x + 1/(x-1) + 1/(x-y))*diff(f,x)        */
  14 /*      - (c_1/x + c_2/(x-1) + c_3/(x-y))* f(x) = 0               */
  15 /*                                                                */
  16 /* and 4 "boundary" conditions: diff(f,x)(a_k)=2*c_k*f(a_k) for   */
  17 /* k=1,2,3, a_k=(0, 1 , y)  and normalization f(1) = 1            */
  18 /*                                                                */
  19 /* The z-vector is z_1=f, z_2=diff(f,x), z_3=c_2 and z_4=c_3      */
  20 /* The guess is chosen to have one node in [0,1],  f(x)=2*x-1     */
  21 /* such that f(1)=1, c_2 and c_3 are chosen to cancel poles in    */
  22 /* the differential equation at 1.0 and y, z_3=1, z_4=1/(2*y-1)   */
  23 /* Ref: http://arxiv.org/pdf/hep-th/0407005                       */
  24 /******************************************************************/
  25
  26
  27
  28
  29 /* Number of differential equations */
  30 ncomp : 3;
  31 /* Orders */
  32 m : [2, 1, 1];
  33
  34
  35 /* Set up parameter array.  */
  36 ipar : makelist(0,k,1,11);
  37
  38 /* Non-linear prob */
  39 ipar[1] : 1;
  40
  41 /* Number of collocation points */
  42 ipar[2] : 3;
  43
  44
  45 /* Initial uniform mesh of 4 subintervals */
  46 ipar[3] : 4;
  47 ipar[8] : 0;
  48
  49 /* Size of fspace, ispace */
  50 ipar[5] : 30000;
  51 ipar[6] : 3000;
  52
  53 /* Medium output */
  54 ipar[7] : 0;
  55 /* Initial approx is provided */
  56 ipar[9] : 1;
  57
  58 /* fixpnt is an array containing all fixed points in the mesh, in
  59 particular "boundary" points, except aleft and aright, ipar[11] its size */
  60 ipar[11] : 1;
  61
  62 /* Tolerances on all components z_1, z_2, z_3, z_4 */
  63 ipar[4] : 4;
  64
  65 /* Tolerance on f and diff(f,x) and on c_2 and c_3 */
  66 ltol : [1, 2, 3, 4];
  67 tol : [1d-5, 1d-5, 1d-5, 1d-5];
  68
  69 fspace : makelist(0d0, k, 1, 30000)$
  70 ispace : makelist(0, k, 1, 3000)$
  71
  72 /* the only one interior "boundary point" */
  73 fixpnt : [1.0d0];
  74
  75
  76 /* Define the differential equations */
  77
  78
  79 /* v is the parameter. To solve the eigenvalues problem one adds the
  80 equations c'_2 = 0, c'_3 = 0, so 4 unknowns constants, and 4 boundary conditions
  81 f -> z[1], f' -> z[2], c_2 -> z[3], c_3 -> z[4], c_1 -> v - z[3] - z[4] */
  82
  83
  84 f(x,z1,z2,z3,z4) :=
  85   [- (1/2.0d0)*(1/x + 1/(x-1.0d0) + 1/(x-y))*z2 +
  86       ((vv - z3 - z4)/x + z3/(x-1.0d0) + z4/(x-y))*z1,
  87   0.0d0, 0.0d0];
  88 df : jacobian(f(x,z1,z2,z3,z4),[z1,z2,z3,z4]);
  89
  90
  91 /* Define the boundary conditions, in order, one at 0, two at 1, one at y */
  92
  93 g(z1,z2,z3,z4) := [z2 - 2.0d0*z1*(vv-z3-z4), z2 - 2.0d0*z1*z3,
  94                 z1 - 1.0d0, z2 - 2.0d0*z1*z4];
  95 dg : jacobian(g(z1,z2,z3,z4), [z1,z2,z3,z4]);
  96
  97
  98 /* Define constant, y is defined here to have clean formulas above */
  99
 100 y: 1.9d0;  /* Can be changed, but effect is not big */
 101 v: 0.00001d0; /* Take v small at first */
 102
 103
 104 /* Locations of side conditions, sorted, depends on y! */
 105 zeta : float([0.0d0, 1.0d0, 1.0d0, y]);
 106 /* Interval ends */
 107 aleft : 0.0d0;
 108 aright : y;
 109
 110
 111 fsub(x,z1,z2,z3,z4) := subst(['x=x,'y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4],
 112                               f(x,z1,z2,z3,z4))$
 113 dfsub(x,z1,z2,z3,z4) := subst(['x=x,'y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4], df)$
 114 gsub(i,z1,z2,z3,z4) := subst(['x=x,'y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4],
 115                               g(z1,z2,z3,z4)[i])$
 116 dgsub(i,z1,z2,z3,z4) := subst(['y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg,i)[1])$
 117
 118
 119
 120 /* eigenstates are characterized by number of nodes in [0,1] and in
 121 [1,y], here guess selects one node (zero) in [0,1] with linear
 122 f(x)=2*x-1 and constant c_2, c_3, so dmval=0. Notice that the z-vector
 123 has mstar = 4 components, while dmval has ncomp = 3 components. */
 124
 125 init_guess(x) := [[2.0d0*x -1.0d0, 2.0d0, 1.0d0, 1/(2.0d0*y-1.0d0)],[0.0d0, 0.0d0, 0.0d0]];
 126
 127
 128
 129 [iflag, fspace, ispace] :
 130   colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
 131                 0, fsub, dfsub, gsub, dgsub, init_guess)$
 132
 133 /* Plot values of solution at x = 0.1 0.2 ... 1.0 ...y */
 134
 135
 136 len: floor(1+y/0.1)$
 137
 138 xx: makelist(0.1*k, k, 0, len)$
 139
 140 zvals: colnew_appsln(xx, 4, fspace, ispace)$
 141
 142 /* Values of the wave function f */
 143
 144 fvals: makelist(zvals[k][1],k, 1, 1+len)$
 145
 146 plot2d([discrete,xx,fvals])$
 147
 148
 149 /* Show eigenvalues, it suffices to use the first value of x. Energy
 150 is 2*(c_2+y*c_3). Show it */
 151
 152 printf(true, "c_2 = ~5,2f    c_3 = ~5,2f   E/2 = ~5,2f ~%",
 153 zvals[1][3], zvals[1][4], zvals[1][3]+y*zvals[1][4])$
 154
 155
 156 /* Now do the continuation by setting ipar[9]=3 and slowly varying v in 40
 157 steps. We store c_2 and c_3 in the array eigenval. */
 158
 159 nsteps:100;
 160 eigenval: makelist([0.0,0.0,0.0],k,0,nsteps)$
 161 eigenval[1][1]: zvals[1][3]$
 162 eigenval[1][2]: zvals[1][4]$
 163 eigenval[1][3]: zvals[1][3]+y*zvals[1][4]$
 164
 165 ipar[9] : 3;
 166 ipar[3] : 1;
 167 for i from 1 thru nsteps do
 168   (v: v + 0.2,
 169    print("v = ", v),
 170    fsub(x,z1,z2,z3,z4) := subst(['x=x,'y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4],
 171                               f(x,z1,z2,z3,z4)),
 172    dfsub(x,z1,z2,z3,z4) := subst(['x=x,'y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4], df),
 173    gsub(i,z1,z2,z3,z4) := subst(['x=x,'y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4],
 174                               g(z1,z2,z3,z4)[i]),
 175    dgsub(i,z1,z2,z3,z4) := subst(['y=y,'vv=v,'z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg,i)[1]),
 176    ipar[3] : ispace[1],
 177    [iflag, fspace, ispace] :
 178    colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
 179                 0, fsub, dfsub, gsub, dgsub, init_guess),
 180    zval: colnew_appsln([0.1d0], 4, fspace, ispace)[1],
 181    eigenval[i+1][1]: zval[3],
 182    eigenval[i+1][2]: zval[4],
 183    eigenval[i+1][3]: zval[3]+y*zval[4])$
 184
 185
 186 /* Plot evolution of eigenvalues with v */
 187
 188 vv: makelist(0.2*k,k,0,nsteps)$
 189 plot2d([[discrete,vv, makelist(eigenval[k][1],k,1,nsteps+1)],
 190   [discrete,vv, makelist(eigenval[k][2],k,1,nsteps+1)],
 191   [discrete,vv, makelist(eigenval[k][3],k,1,nsteps+1)]])$