share/colnew/prob2.mac

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