share/colnew/prob1.mac

   1 /* Problem 1 for colnew */
   2
   3 /* One differential equation of order 4 */
   4 m : [4];
   5
   6 /* Location of boundary conditions */
   7 zeta : float([1,1,2,2]);
   8
   9 /* Set up parameter array.  Use defaults for all except initial mesh
  10 size, number of tolerances and sizes of work arrays */
  11 ipar : makelist(0,k,1,11);
  12 ipar[3] : 1;
  13 ipar[4] : 2;
  14 ipar[5] : 2000;
  15 ipar[6] : 200;
  16
  17 /* Two error tolerances (on u and its second derivative */
  18 ltol : [1, 3];
  19 tol : [1d-7, 1d-7];
  20
  21 fspace : makelist(0d0, k, 1, 2000)$
  22 ispace : makelist(0, k, 1, 200)$
  23 fixpnt : [];
  24
  25 /* Define the equations */
  26 fsub(x, z1, z2, z3, z4) := [(1-6*x^2*z4-6*x*z3)/x^3];
  27 df : jacobian(fsub(x,z1, z2, z3, z4),[z1,z2,z3,z4]);
  28 dfsub(x, z1, z2, z3, z4) := subst(['x=x,'z1=z1,'z2=z2,'z3=z3,'z4=z4], df);
  29 g(z1, z2, z3, z4) := [z1, z3, z1, z3];
  30 gsub(i, z1, z2, z3, z4) := subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], g(z1, z2, z3, z4)[i]);
  31 dg:jacobian(g(z1, z2, z3, z4), [z1,z2,z3,z4]);
  32 dgsub(i, z1, z2, z3, z4) := subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg, i)[1]);
  33
  34 exact(x) := [.25*(10.*log(2.)-3.)*(1.-x) + .5*(1./x+(3.+x)*log(x)-x),
  35              -.25*(10.*log(2.)-3.) + .5*(-1./x/x+log(x)+(3.+x)/x-1.),
  36              .5*(2./x**3+1./x-3./x/x),
  37              .5*(-6./x**4-1./x/x+6./x**3)];
  38 [iflag, fspace, ispace] :
  39   colnew_expert(1, m, 1d0, 2d0, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
  40                 0, fsub, dfsub, gsub, dgsub, dummy)$
  41
  42 /* Calculate the error at 101 points using the known exact solution */
  43
  44 block([x : 1,
  45        err : makelist(0d0, k, 1, 4),
  46        j],
  47   for j : 1 thru 101 do
  48     block([],
  49       zval : colnew_appsln([x], 4, fspace, ispace)[1],
  50       u : float(exact(x)),
  51       err : map(lambda([a,b], max(a,b)), err, abs(u-zval)),
  52       x : x + 0.01),
  53   print("The exact errors are:"),
  54   printf(true, "   ~{ ~11,4e~}~%", err));