share/colnew/prob3.mac

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