share/colnew/prob2.mac

   1 /* Problem 2 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 /* Define constants */
  11
  12 gamma : 1.1;
  13 eps : 0.01;
  14 dmu : eps;
  15 eps4mu : eps^4/dmu;
  16 xt : sqrt(2*(gamma-1)/gamma);
  17
  18 /* Number of differential equations */
  19 ncomp : 2;
  20 /* Orders */
  21 m : [2, 2];
  22
  23 /* Interval ends */
  24 aleft : 0.0;
  25 aright : 1.0;
  26
  27 /* Locations of side conditions */
  28 zeta : float([0, 0, 1, 1]);
  29
  30 /* Set up parameter array.  */
  31 ipar : makelist(0,k,1,11);
  32 /* Non-linear prob */
  33 ipar[1] : 1;
  34 /* 4 collocation points per subinterval */
  35 ipar[2] : 4;
  36 /* Initial uniform mesh of 10 subintervals */
  37 ipar[3] : 10;
  38 ipar[8] : 0;
  39 /* Size of fspace, ispace */
  40 ipar[5] : 40000;
  41 ipar[6] : 2500;
  42 /* Full output */
  43 ipar[7] : -1;
  44 /* Initial approx is provided */
  45 ipar[9] : 1;
  46 /* Regular problem */
  47 ipar[10] : 0;
  48 /* No fixed points in mesh */
  49 ipar[11] : 0;
  50 /* Tolerances on all components */
  51 ipar[4] : 4;
  52
  53 /* Two error tolerances (on u and its second derivative */
  54 ltol : [1, 2, 3, 4];
  55 tol : [1d-5, 1d-5, 1d-5, 1d-5];
  56
  57 fspace : makelist(0d0, k, 1, 40000)$
  58 ispace : makelist(0, k, 1, 2500)$
  59 fixpnt : [];
  60
  61 /* Define the equations */
  62 fsub(x, z1, z2, z3, z4) := [z1/x/x - z2/x + (z1-z3*(1-z1/x) - gamma*x*(1-x*x/2))/eps4mu,
  63                z3/x/x - z4/x + z1*(1-z1/2/x)/dmu];
  64 df : jacobian(fsub(x,z1, z2, z3, z4),[z1,z2,z3,z4]);
  65 dfsub(x, z1, z2, z3, z4) :=
  66   subst(['x=x,'z1=z1,'z2=z2,'z3=z3,'z4=z4], df);
  67 g(z1, z2, z3, z4) := [z1, z3, z1, z4 - 0.3*z3 + .7];
  68 gsub(i, z1, z2, z3, z4) :=
  69     subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], g(z1, z2, z3, z4)[i]);
  70 dg:jacobian(g(z1, z2, z3, z4), [z1,z2,z3,z4]);
  71 dgsub(i, z1, z2, z3, z4) := subst(['z1=z1,'z2=z2,'z3=z3,'z4=z4], row(dg, i)[1]);
  72
  73 solutn(x) :=
  74   block([cons : gamma*x*(1-0.5*x*x),
  75          dcons : gamma*(1-1.5*x*x),
  76          d2cons : -3*gamma*x],
  77     if is(x > xt) then
  78       [[0, 0, -cons, -dcons],
  79        [0, -d2cons]]
  80     else
  81       [[2*x, 2, -2*x + cons, -2 + dcons],
  82        [0, d2cons]]);
  83
  84 [iflag, fspace, ispace] :
  85   colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace,
  86                 0, fsub, dfsub, gsub, dgsub, solutn);
  87
  88 /* Print values of solution at x = 0,.05,...,1 */
  89 x : 0;
  90 for j : 1 thru 21 do
  91   block([],
  92     zval : colnew_appsln([x], 4, fspace, ispace)[1],
  93     printf(true, "~5,2f  ~{~15,5e~}~%", x, zval),
  94     x : x + 0.05);