share/colnew/prob5.mac

   1 /* This example (Ascher et al, 1995, Example 9.2) solves a numerically
   2    difficult boundary value problem using continuation.
   3
   4 Reference:
   5   Uri M. Ascher, Robert M. M. Mattheij and Robert D. Russell.
   6   Numerical Solution of Boundary Value Problems for Ordinary
   7   Differential Equations. Classics in Applied Mathematics 13, SIAM, (1995),
   8   ISBN 978-0-89871-354-1
   9
  10 The linear differential equation is
  11
  12     e u''(x) + x u'(x) = -e %pi^2* cos(%pi x) - (%pi x) sin(%pi x)
  13
  14 with -1 < x < 1 and 0 < e <= 1
  15
  16 and boundary conditions u(-1)=-2 and u(1)=0
  17
  18 The exact solution is
  19     u(x) = cos(%pi x) + erf(x/sqrt(2 e))/erf(1/sqrt(2 e))
  20
  21 When e is small the solution has a rapid transition near x=0
  22 and is difficult to solve numerically.  COLNEW is able to solve the
  23 problem for directly for e=1.0e-6, but here we will use
  24 continuation to solve it successively for e=[1e-2,1e-3,1e-4,1e-5,1e-6].
  25
  26 The mesh and solution from the previous step is used as the initial
  27 estimates for the next step.
  28
  29 */
  30
  31 load("colnew")$
  32
  33 kill(e,x,z1,z2)$
  34
  35 /* Exact solution */
  36 exact(x):=cos(%pi*x)+erf(x/sqrt(2*e))/erf(1/sqrt(2*e))$
  37
  38 /* Define the equations.  Do this before e is defined */
  39 f: [-(x/e)*z2 - %pi^2*cos(%pi*x) - (%pi*x/e)*sin(%pi*x)];
  40 define(fsub(x,z1,z2),f);
  41 df: jacobian(f,[z1,z2]);
  42 define(dfsub(x,z1,z2),df);
  43
  44 /* Build the functions gsub and dgsub
  45    Use define and buildq to remove dependence on g and dg */
  46 g: [z1+2,z1];
  47 define(gsub(i,z1,z2),buildq([g],g[i]));
  48 dg: jacobian(g,[z1,z2]);
  49 define(
  50   dgsub(i,z1,z2),
  51   buildq([val:makelist(dg[i],i,1,length(dg))],block([dg:val],dg[i])));
  52
  53 /* Define constant epsilon */
  54 e : 0.01$
  55
  56 /* Number of differential equations */
  57 ncomp : 1$
  58 /* Orders */
  59 m : [2]$
  60
  61 /* Interval ends */
  62 aleft  : -1.0$
  63 aright :  1.0$
  64
  65 /* Locations of side conditions */
  66 zeta : float([-1, 1])$
  67 /* Set up parameter array.  */
  68 ipar : makelist(0,k,1,11)$
  69 /* linear prob */
  70 ipar[1] : 0$
  71 /* 5 collocation points per subinterval */
  72 ipar[2] : 5$
  73 /* Initial uniform mesh of 1 subintervals */
  74 ipar[3] : 1$
  75 ipar[8] : 0$
  76 /* Size of fspace, ispace */
  77 ipar[5] : 10000$
  78 ipar[6] :  1000$
  79 /* Minimal output */
  80 ipar[7]:0$
  81 /* No initial guess is provided */
  82 ipar[9] : 0$
  83 /* Regular problem */
  84 ipar[10] : 0$
  85 /* No fixed points in mesh */
  86 ipar[11] : 0$
  87 /* Tolerances on two components */
  88 ipar[4] : 2$
  89
  90 /* Two error tolerances (on u and its derivative)
  91    Relatively large tolerances to keep the example small */
  92 ltol : [1, 2]$
  93 tol : [1e-4, 1e-4]$
  94
  95 fspace : makelist(0e0, k, 1, ipar[5])$
  96 ispace : makelist(0, k, 1, ipar[6])$
  97 fixpnt : []$
  98
  99 /* First run with default initial guess.
 100    Returns iflag. 1 = success */
 101 ([iflag, fspace, ispace] :
 102 colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol,
 103 fixpnt, ispace, fspace,
 104 0, fsub, dfsub, gsub, dgsub, dummy), iflag);
 105
 106 /* iflag=1 on success */
 107 if (iflag#1) then error("On return from colnew_expert: iflag # 1");
 108
 109 /* Function to generate equally spaced list of values */
 110 xlist(xmin,xmax,n):=block([dx:(xmax-xmin)/n],makelist(i,i,0,n)*dx+xmin)$
 111
 112 /* x values for solution.  Cluster around x=0 */
 113 X: xlist(aleft,aright,500)$
 114
 115 /* Generate solution values for z1=u(x) */
 116 ans:colnew_appsln(X,2,fspace,ispace)$
 117 z:maplist(first,ans)$
 118 Z:[z]$
 119
 120 /* Compare with exact solution and report */
 121 y:float(map(exact,X))$
 122 maxerror:apply(max,abs(y-z));
 123 printf(true," e: ~8,3e  iflag ~3d  Mesh size ~3d  max error ~8,3e~%",
 124   e,iflag,ispace[1],maxerror);
 125
 126 /* Now use continuation to solve for progressively smaller e
 127    Use previous solution as initial guess and every second point
 128    from previous mesh as initial mesh */
 129 ipar[9] : 3$
 130
 131 /* Run COLNEW using continuation for new value of e
 132    Set new mesh size ipar[3] from previous size ispace[1]
 133    Push list of values of z1=u(x) on to list Z */
 134 run_it(e_):=block(
 135   e:e_,
 136   ipar[3]:ispace[1],
 137   [iflag, fspace, ispace]:
 138      colnew_expert(ncomp,m,aleft,aright,zeta,ipar,ltol,tol,fixpnt,
 139      ispace,fspace,0,fsub,dfsub,gsub,dgsub,dummy),
 140   if (iflag#1) then error("On return from colnew_expert: iflag # 1"),
 141   ans:colnew_appsln(X,2,fspace,ispace),
 142   z:maplist(first,ans),
 143   push(z,Z),
 144   y:float(map(exact,X)),
 145   maxerror:apply(max,abs(y-z)),
 146   printf(true," e: ~8,3e  iflag ~3d  Mesh size ~3d  max error ~8,3e~%",
 147     e,iflag,ispace[1],maxerror),
 148   iflag
 149 );
 150
 151 for e_ in [1e-3,1e-4,1e-5,1e-6] do run_it(e_);
 152
 153 /* Z is list of solutions z1 = u(x).  Restore order. */
 154 Z:reverse(Z)$
 155
 156 /* Plot z1=u(x) for each value of e*/
 157 /*
 158 plot2d([
 159   [discrete,X,Z[1]], [discrete,X,Z[2]], [discrete,X,Z[3]],
 160   [discrete,X,Z[4]], [discrete,X,Z[5]]],
 161   [legend,"e=1e-2","e=1e-3","e=1e-4","e=1e-5","e=1e-6"],
 162   [xlabel,"x"],[ylabel,"u(x)"],
 163   [png_file,"./colnew-ex5.png"]);
 164 */