tests/rtestode_zp.mac

   1 /* Tests from:
   2     F Postel and P Zimmermann, A Review of the ODE solvers of
   3     AXIOM, DERIVE, MACSYMA, MATHEMATICA, MUPAD and REDUCE
   4     Proceedings of the 5th Rhine Workshop on Computer Algebra
   5     April 1-3, 1996, Saint-Louis, France
   6     http://www.loria.fr/~zimmerma/ComputerAlgebra/
   7 */
   8
   9 kill(all);
  10 done;
  11
  12 /* equation 1 - linear polynomial */
  13 eqn: (x^4-x^3)*diff(u(x),x) + 2*x^4*u(x) = x^3/3 + c;
  14 (x^4-x^3)*'diff(u(x),x,1)+2*x^4*u(x) = x^3/3+c;
  15 ans:ode2(eqn,u(x),x);
  16 u(x) = ((2*x^3-3*x^2+6*c)*%e^(2*x)/(12*x^2)+%c)*%e^-(2*(log(x-1)+x));
  17 factor(ans);
  18 u(x)=%e^-(2*x)*(2*x^3*%e^(2*x)-3*x^2*%e^(2*x)+6*c*%e^(2*x)+12*%c*x^2)/(12*(x-1)^2*x^2);
  19 is(ratsimp(ev(eqn,ans,diff)));
  20 true;
  21
  22 /* equation (2) {firstord2} */
  23 eqn: -1/2*'diff(y,x) + y = sin(x);
  24 y-'diff(y,x,1)/2 = sin(x);
  25 ans: ode2(eqn,y,x);
  26 y = %e^(2*x)*(%c-2*%e^-(2*x)*(-2*sin(x)-cos(x))/5);
  27 ans: expand(ans);
  28 y = 4*sin(x)/5+2*cos(x)/5+%c*%e^(2*x);
  29 is(ratsimp(ev(eqn,ans,diff)));
  30 true;
  31
  32 /* equation 3 - linear change of variables */
  33 eqn: 'diff(y,x,2)*(a*x+b)^2+4*'diff(y,x)*(a*x+b)*a+2*y*a^2=0;
  34 (a*x+b)^2*'diff(y,x,2)+4*a*(a*x+b)*'diff(y,x,1)+2*a^2*y = 0;
  35 ans:ode2(eqn,y,x);
  36 y = %k1*x/(a*x+b)^2+%k2/(a*x+b)^2;
  37 is(ratsimp(ev(eqn,%,diff)));
  38 true;
  39
  40 /* equation 4 - linear polynomial (adjoint equation)
  41    see zwillinger, p151
  42    this crashed maxima-5.5
  43    */
  44 eqn: (x^2-x)*'diff(y,x,2)+(2*x^2+4*x-3)*'diff(y,x)+8*x*y=1;
  45 (x^2-x)*'diff(y,x,2)+(2*x^2+4*x-3)*'diff(y,x,1)+8*x*y = 1;
  46 ans:ode2(eqn,y,x);
  47 false;
  48
  49 /* equation 5 - linear polynomial */
  50 eqn: (x^2-x)*'diff(y,x,2)+(1-2*x^2)*'diff(y,x)+(4*x-2)*y=0;
  51 (x^2-x)*'diff(y,x,2)+(1-2*x^2)*'diff(y,x,1)+(4*x-2)*y = 0;
  52 ans:ode2(eqn,y,x);
  53 false;
  54
  55 /* equation 6 - dependent variable missing */
  56 eqn:'diff(y,x,2)+2*x*'diff(y,x)=2*x;
  57 'diff(y,x,2)+2*x*'diff(y,x,1) = 2*x;
  58 ans: ode2(eqn,y,x);
  59 y = sqrt(%pi)*%k1*erf(x)/2+x+%k2;
  60 is(ratsimp(ev(eqn,ans,diff)));
  61 true;
  62
  63 /* equation 7 - liouvillian solution */
  64 eqn: (x^3/2-x^2)*'diff(y,x,2)+(s2*x^2-3*x+1)*'diff(y,x)+(x-1)*y=0;
  65 (x^3/2-x^2)*'diff(y,x,2)+(s2*x^2-3*x+1)*'diff(y,x,1)+(x-1)*y = 0;
  66 ode2(eqn,y,x);
  67 false;
  68
  69 /* equation 8 - reduction of order */
  70 eqn: 'diff(y,x,2)-2*x*'diff(y,x)+2*y=3;
  71 'diff(y,x,2)-2*x*'diff(y,x)+2*y = 3;
  72 ode2(eqn,y,x);
  73 false;
  74
  75 /* equation 9 - integrating factors */
  76 eqn: sqrt(x)*'diff(y,x,2)+2*x*'diff(y,x)+3*y=0;
  77 sqrt(x)*'diff(y,x,2)+2*x*'diff(y,x)+3*y=0;
  78 ode2(eqn,y,x);
  79 false;
  80
  81 /* equation 18 - bernoulli equation */
  82 eqn: 'diff(y,x) + y = y^3*sin(x);
  83 'diff(y,x,1)+y = sin(x)*y^3;
  84 ans: ode2(eqn,y,x);
  85 y = %e^-x/sqrt(%c-2*%e^-(2*x)*(-2*sin(x)-cos(x))/5);
  86 is(ratsimp(ev(eqn,ans,diff)));
  87 true;
  88
  89 /* equation (20) {clairaut} */
  90 eqn: (x^2-1)*'diff(y,x)^2 - 2*x*y*'diff(y,x) + y^2 -1 = 0;
  91 (x^2-1)*'diff(y,x)^2 - 2*x*y*'diff(y,x) + y^2 -1 = 0;
  92 ode2(eqn,y,x);
  93 false;