share/colnew/prob1.mac

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