## maxima-commits

 [Maxima-commits] CVS: maxima/share/colnew prob1.mac, NONE, 1.1 prob2.mac, NONE, 1.1 From: Raymond Toy - 2010-05-27 18:33:40 ```Update of /cvsroot/maxima/maxima/share/colnew In directory sfp-cvsdas-4.v30.ch3.sourceforge.com:/tmp/cvs-serv22057/share/colnew Added Files: prob1.mac prob2.mac Log Message: prob1.mac: prob2.mac: o Simple demo of colnew. --- NEW FILE: prob1.mac --- /* Problem 1 for colnew */ /* One differential equation of order 4 */ m : [4]; /* Location of boundary conditions */ zeta : float([1,1,2,2]); /* Set up parameter array. Use defaults for all except initial mesh size, number of tolerances and sizes of work arrays */ ipar : makelist(0,k,1,11); ipar[3] : 1; ipar[4] : 2; ipar[5] : 2000; ipar[6] : 200; /* Two error tolerances (on u and its second derivative */ ltol : [1, 3]; tol : [1d-7, 1d-7]; fspace : makelist(0d0, k, 1, 2000)\$ ispace : makelist(0, k, 1, 200)\$ fixpnt : []; /* Define the equations */ fsub(x, z) := [(1-6*x^2*z[4]-6*x*z[3])/x^3]; df : jacobian(fsub(x,z),[z[1],z[2],z[3],z[4]]); dfsub(x, z) := ''df; g(z) := [z[1], z[3], z[1], z[3]]; gsub(i, z) := g(z)[i]; dg:jacobian(g(z), [z[1],z[2],z[3],z[4]]); dgsub(i, z) := row(dg, i)[1]; exact(x) := [.25*(10.*log(2.)-3.)*(1.-x) + .5*(1./x+(3.+x)*log(x)-x), -.25*(10.*log(2.)-3.) + .5*(-1./x/x+log(x)+(3.+x)/x-1.), .5*(2./x**3+1./x-3./x/x), .5*(-6./x**4-1./x/x+6./x**3)]; [iflag, fspace, ispace] : colnew_expert(1, m, 1d0, 2d0, zeta, ipar, ltol, tol, fixpnt, ispace, fspace, 0, fsub, dfsub, gsub, dgsub, dummy)\$ /* Calculate the error at 101 points using the known exact solution */ block([x : 1, err : makelist(0d0, k, 1, 4), j], for j : 1 thru 101 do block([], zval : colnew_appsln([x], 4, fspace, ispace)[1], u : float(exact(x)), err : map(lambda([a,b], max(a,b)), err, abs(u-zval)), x : x + 0.01), print("The exact errors are:"), printf(true, " ~{ ~11,4e~}~%", err)); --- NEW FILE: prob2.mac --- /* Problem 2 for colnew */ /* Define constants */ gamma : 1.1; eps : 0.01; dmu : eps; eps4mu : eps^4/dmu; xt : sqrt(2*(gamma-1)/gamma); /* Number of differential equations */ ncomp : 2; /* Orders */ m : [2, 2]; /* Interval ends */ aleft : 0.0; aright : 1.0; /* Locations of side conditions */ zeta : float([0, 0, 1, 1]); /* Set up parameter array. */ ipar : makelist(0,k,1,11); /* Non-linear prob */ ipar[1] : 1; /* 4 collocation points per subinterval */ ipar[2] : 4; /* Initial uniform mesh of 10 subintervals */ ipar[3] : 10; ipar[8] : 0; /* Size of fspace, ispace */ ipar[5] : 40000; ipar[6] : 2500; /* Full output */ ipar[7] : -1; /* Initial approx is provided */ ipar[9] : 1; /* Regular problem */ ipar[10] : 0; /* No fixed points in mesh */ ipar[11] : 0; /* Tolerances on all components */ ipar[4] : 4; /* Two error tolerances (on u and its second derivative */ ltol : [1, 2, 3, 4]; tol : [1d-5, 1d-5, 1d-5, 1d-5]; fspace : makelist(0d0, k, 1, 40000)\$ ispace : makelist(0, k, 1, 2500)\$ fixpnt : []; /* Define the equations */ fsub(x, z) := [z[1]/x/x - z[2]/x + (z[1]-z[3]*(1-z[1]/x) - gamma*x*(1-x*x/2))/eps4mu, z[3]/x/x - z[4]/x + z[1]*(1-z[1]/2/x)/dmu]; df : jacobian(fsub(x,z),[z[1],z[2],z[3],z[4]]); dfsub(x, z) := ''df; g(z) := [z[1], z[3], z[1], z[4] - 0.3*z[3] + .7]; gsub(i, z) := g(z)[i]; dg:jacobian(g(z), [z[1],z[2],z[3],z[4]]); dgsub(i, z) := row(dg, i)[1]; solutn(x) := block([cons : gamma*x*(1-0.5*x*x), dcons : gamma*(1-1.5*x*x), d2cons : -3*gamma*x], if is(x > xt) then [[0, 0, -cons, -dcons], [0, -d2cons]] else [[2*x, 2, -2*x + cons, -2 + dcons], [0, d2cons]]); [iflag, fspace, ispace] : colnew_expert(ncomp, m, aleft, aright, zeta, ipar, ltol, tol, fixpnt, ispace, fspace, 0, fsub, dfsub, gsub, dgsub, solutn); /* Print values of solution at x = 0,.05,...,1 */ x : 0; for j : 1 thru 21 do block([], zval : colnew_appsln([x], 4, fspace, ispace)[1], printf(true, "~5,2f ~{~15,5e~}~%", x, zval), x : x + 0.05); ```