tests/rtestode.mac

   1 /* ODE tests */
   2
   3 kill(all);
   4 done;
   5
   6 /* Bug #2796 ode2 with n declared constant
   7  *
   8  * Can't load ode2.mac if any formal function args are
   9  * declared constant
  10  */
  11 (declare(n,constant), load(ode2), remove(n,constant));
  12 done;
  13
  14 /* Trivial ode - bug 866510 */
  15 ode2('diff(y,x),y,x);
  16 y=%c;
  17
  18 /* Examples from "The Maxima Book" */
  19
  20 ode2(x^2*'diff(y,x)+3*x*y=sin(x)/x, y, x);
  21 y = (%c-cos(x))/x^3;
  22 ic1(%, x=1, y=1);
  23 y = -((cos(x)-cos(1)-1)/x^3);
  24 method;
  25 linear;
  26
  27 soln:ode2('diff(y,x,2) + y = 4*x, y, x);
  28 y = %k1*sin(x) + %k2*cos(x) + 4*x;
  29 method;
  30 variationofparameters;
  31 ic2(soln, x=0, y=1, 'diff(y,x)=3);
  32 y = -sin(x)+cos(x)+4*x;
  33 bc2(soln, x=0, y=3, x=2, y=1);
  34 y = -((3*cos(2)+7)*sin(x)/sin(2)) + 3*cos(x) + 4*x;
  35
  36 ode2((3*x^2+4*x+2)=(2*y-1)*'diff(y,x), y, x);
  37 y^2-y = x^3+2*x^2+2*x+%c;
  38 method;
  39 separable;
  40
  41 ode2(x^2*cos(x*y)*'diff(y,x) + (sin(x*y)+x*y*(cos(x*y)))=0, y, x);
  42 x*sin(x*y)=%c;
  43 method;
  44 exact;
  45
  46 ode2( (2*x*y-exp(-2*y))*'diff(y,x)+y=0, y, x);
  47 x*exp(2*y) - log(y) = %c;
  48 method;
  49 exact;
  50 intfactor;
  51 exp(2*y)/y;
  52
  53 ode2( 'diff(y,x)=(y/x)^2+2*(y/x), y, x);
  54 -((x*y+x^2)/y) = %c;
  55 method;
  56 exact;
  57
  58 ode2( 'diff(y,x)+(2/x)*y=(1/x^2)*y^3, y, x);
  59 y = 1/(sqrt( 2/(5*x^5) + %c)*x^2);
  60 method;
  61 bernoulli;
  62 odeindex;
  63 3;
  64
  65 ode2( 'diff(y,x,2)-3*'diff(y,x)+2*y=0, y, x);
  66 y = %k1*exp(2*x) + %k2*exp(x);
  67 method;
  68 constcoeff;
  69
  70 ode2( 'diff(y,x,2)-4*'diff(y,x)+4*y=0, y, x);
  71 y = (%k2*x + %k1)*exp(2*x);
  72 method;
  73 constcoeff;
  74
  75 ode2(x^2*'diff(y,x,2)+x*'diff(y,x)-y=0, y, x);
  76 y=%k2*x-%k1/(2*x);
  77 method;
  78 exact;
  79
  80 ode2( x^2*'diff(y,x,2)+4*x*'diff(y,x)+2*y=0, y, x);
  81 y=%k1/x+%k2/x^2;
  82 method;
  83 exact; /*euler*/
  84
  85 ode2( x^2*'diff(y,x,2)+5*x*'diff(y,x)+4*y=0, y, x);
  86 y=(%k2*log(x)+%k1)/x^2;
  87 method;
  88 euler;
  89
  90 ode2( x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-1/4)*y=0, y, x);
  91 y=(%k1*sin(x)+%k2*cos(x))/sqrt(x);
  92 method;
  93 bessel;
  94
  95 ode2( x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-4)*y=0, y, x);
  96 y=%k1*bessel_j(2,x)+%k2*bessel_y(2,x);
  97 method;
  98 bessel;
  99
 100 ode2( (x-1)^2*'diff(y,x,2)+(x-1)*'diff(y,x)+((x-1)^2-4)*y=0, y, x);
 101 y=%k1*bessel_j(2,x-1)+%k2*bessel_y(2,x-1);
 102 method;
 103 bessel;
 104
 105 ode2( x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-1/9)*y=0, y, x);
 106 y=bessel_j(-1/3,x)*%k2+bessel_j(1/3,x)*%k1;
 107 method;
 108 bessel;
 109
 110 /* Bug report 2876387: asks if obvious non-integers are integers */
 111 (declare(n,integer),ode2(x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-n^2)*y=0,y,x));
 112 y = %k2*bessel_y(n,x)+%k1*bessel_j(n,x);
 113 (remove(n,integer),method);
 114 bessel;
 115
 116 (declare(v,noninteger),ode2(x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-v^2)*y=0,y,x));
 117 y = %k1*bessel_j(v,x)+%k2*bessel_j(-v,x);
 118 (remove(v,noninteger),method);
 119 bessel;
 120
 121 ode2(x^2*'diff(y,x,2)+x*'diff(y,x)+(x^2-3)*y=0,y,x);
 122 y = %k1*bessel_j(sqrt(3),x)+%k2*bessel_j(-sqrt(3),x);
 123 method;
 124 bessel;
 125
 126 ode2( 'diff(y,x,2)+2*'diff(y,x)+y=exp(x), y, x);
 127 y=exp(x)/4+(%k2*x+%k1)*exp(-x);
 128 method;
 129 variationofparameters;
 130 yp;
 131 exp(x)/4;
 132
 133 ode2( x*'diff(y,x,2)+('diff(y,x))^2=0, y, x);
 134 /* y='integrate(1/(log(x)+%k1),x)+%k2;
 135    Because of adding more integrals for the power function we get a result
 136    12/2008 */
 137 y=%k2-expintegral_e(1,-log(x)-%k1)*%e^-%k1;
 138 method;
 139 freeofy;
 140
 141 ode2( y*'diff(y,x,2)+('diff(y,x))^2=0, y, x);
 142 y^2/(2*%k1)=x+%k2;
 143 method;
 144 freeofx;
 145
 146 eq: 'diff(y,x,2)+x*'diff(y,x)+exp(-x^2)*y=0;
 147 'diff(y,x,2)+x*'diff(y,x,1)+%e^-x^2*y = 0;
 148 ans:ode2(eq,y,x);
 149 y = %k1*sin((1/2) * sqrt(2)*sqrt(%pi)*erf(x/sqrt(2)))+%k2*cos((1/2) * sqrt(2)*sqrt(%pi)*erf(x/sqrt(2)));
 150 is(ratsimp(ev(eq,ans,diff)));
 151 true;
 152 method;
 153 xformtoconstcoeff;
 154
 155 eq:x*'diff(y,x,2)+(x^2-1)*'diff(y,x,1)+x^3*y=0;
 156 x*'diff(y,x,2)+(x^2-1)*'diff(y,x,1)+x^3*y=0;
 157 ans:ode2(eq,y,x);
 158 y=%e^-(x^2/4)*(%k1*sin(sqrt(3)*x^2/4)+%k2*cos(sqrt(3)*x^2/4));
 159 is(ratsimp(ev(eq,ans,diff)));
 160 true;
 161 method;
 162 xformtoconstcoeff;
 163
 164 /* Tests of desolve */
 165
 166 eqn1:'diff(f(x),x) = sin(x)+'diff(g(x),x);
 167 'diff(f(x),x,1) = 'diff(g(x),x,1)+sin(x);
 168 eqn2:'diff(g(x),x,2) = 'diff(f(x),x)-cos(x);
 169 'diff(g(x),x,2) = 'diff(f(x),x,1)-cos(x);
 170 desolve([eqn1,eqn2],[f(x),g(x)]);
 171 [f(x)=%e^x*(at('diff(g(x),x,1),x = 0))-at('diff(g(x),x,1),x = 0)+f(0),g(x)=%e^x*(at('diff(g(x),x,1),x=0))-at('diff(g(x),x,1),x = 0)+cos(x)+g(0)-1];
 172 atvalue('diff(g(x),x),x = 0,a);
 173 a;
 174 atvalue(f(x),x = 0,1);
 175 1;
 176 desolve([eqn1,eqn2],[f(x),g(x)]);
 177 [f(x) = a*%e^x-a+1,g(x) = cos(x)+a*%e^x-a+g(0)-1];
 178 remove(f,atvalue,g,atvalue);
 179 done;
 180
 181 atvalue('diff(g(x),x),x = 0,a);
 182 a;
 183 atvalue(f(x),x = 0,1);
 184 1;
 185 desolve([eqn1,eqn2],[f(x),g(x)]);
 186 [f(x) = a*%e^x-a+1,g(x) = cos(x)+a*%e^x-a+g(0)-1];
 187
 188 eqn3: 'diff(f(x),x,2)+f(x)=2*x;
 189 'diff(f(x),x,2)+f(x)=2*x;
 190 desolve(eqn3,f(x));
 191 ''(f(x) = sin(x)*(at('diff(f(x),x,1),x = 0)-2)+f(0)*cos(x)+2*x);
 192
 193 /* BUG 3899 */
 194 errcatch(desolve('diff(f(x),x),x));
 195 [];
 196
 197 /* Examples mentioned in bug report [ 1063454 ] bug in ode2
 198  * First one was reported to fail in CMUCL with "run out of heap" message.
 199  * Others were reported to be OK. Put them all here for good measure.
 200  */
 201
 202 (ode2 ('diff(y, t, 2) + 'diff(y, t) + y - sin(t), y, t),
 203  rhs(%%), ratsimp (diff(%%, t, 2) + diff(%%, t) + %% - sin(t)));
 204 0;
 205
 206 (ode2 ('diff(y, t, 2) + 'diff(y, t) + 2*y - sin(t), y, t),
 207  rhs(%%), ratsimp (diff(%%, t, 2) + diff(%%, t) + 2*%% - sin(t)));
 208 0;
 209
 210 (ode2 ('diff(y, t, 2) + 'diff(y, t) + y - exp(%i*t), y, t),
 211  rhs(%%), ratsimp (diff(%%, t, 2) + diff(%%, t) + %% - exp(%i*t)));
 212 0;
 213
 214 /* bug report 1063454 claims "maxima gets stuck" on the following */
 215 (integrate (my_integrand : exp(t/2) * sin(t) * sin(sqrt(3) * t/2), t),
 216  ratsimp (exponentialize (diff (%%, t) - my_integrand)));
 217 0;
 218
 219 /* Examples to show that ic2 works as expected after revision 1.5 of ode.mac
 220  */
 221
 222 'diff(y,x,2)+y*('diff(y,x,1))^3 = 0;
 223 'diff(y,x,2)+y*('diff(y,x,1))^3 = 0;
 224
 225 soln:ode2(%,y,x);
 226 (y^3+6*%k1*y)/6 = x+%k2;
 227
 228 ratsimp(ic2(soln, x=0, y=0, 'diff(y,x,1)=2));
 229 (y^3+3*y)/6=x;
 230
 231 ratsimp(ic2(soln, x=0, y=0, 'diff(y,x,1)=1));
 232 (y^3+6*y)/6 = x$
 233
 234 /* This is the example of the bug report
 235  * ID:2881021 - ic2 and bc2 may return incorrect results (solution suggeste)
 236  */
 237
 238 ratsimp(ic2(soln, x=0, y=1, 'diff(y,x,1)=2));
 239 y^3/6 = (6*x+1)/6;
 240
 241 /* These examples show that ic2 works for a list of equation and nested
 242  * lists of equation.
 243  */
 244
 245 ratsimp(ic2([soln, soln, soln], x=0, y=0, 'diff(y,x,1)=2));
 246 [(y^3+3*y)/6=x, (y^3+3*y)/6=x, (y^3+3*y)/6=x];
 247
 248 ratsimp(ic2([soln, [soln, soln]], x=0, y=0, 'diff(y,x,1)=2));
 249 [(y^3+3*y)/6=x, [(y^3+3*y)/6=x, (y^3+3*y)/6=x]];
 250
 251 /* Bug report ID: 1839088 - ic2 fails with division by 0
 252  * Maxima no longer gives an error, but does not find the solution.
 253  */
 254
 255 ode2(y*'diff(y,x,2)=a, y, x);
 256 [sqrt(%pi)*%i*%e^-%k1*erf(%i*sqrt(a*log(y)+%k1*a)/sqrt(a))/(sqrt(2)*sqrt(a))
 257    = x+%k2,
 258  -sqrt(%pi)*%i*%e^-%k1*erf(%i*sqrt(a*log(y)+%k1*a)/sqrt(a))/(sqrt(2)*sqrt(a))
 259    = x+%k2];
 260
 261 ic2(%,x=0,y=b,diff(y,x)=0);
 262 [];
 263
 264 /* Bug report ID: 2997443 - ic2 fails
 265  * Maxima no longer gives an error, but does not find the solution.
 266  * The solution of ic2 could be: y=1/20*(sqrt(160*x+1)-1)
 267  */
 268
 269 ode2('diff(x,t,2)+5*'diff(x,t)^3, x, t);
 270 [x = %k2-2*1/sqrt(1/(t+%k1))/sqrt(10),x = 2*1/sqrt(1/(t+%k1))/sqrt(10)+%k2];
 271
 272 ic2(%,t=0,x=0,'diff(x,t)=4);
 273 [];
 274
 275 /* Bug report ID: 1789213 - ic1 for solution containing indefinite integral
 276  * More general implementation of ic1, which handles a noun form of an
 277  * integral correctly. The result simplifies correctly, if we define
 278  * a function and reevaluate the result.
 279  */
 280 sol: ode2(kappa(p) = -'diff(V, p) / V, V, p);
 281 V = %c*%e^-'integrate(kappa(p),p);
 282
 283 ic1(sol, p = p0, V = V0);
 284 V = V0*%e^('at('integrate(kappa(p),p),[p = p0,V = V0])
 285           -'integrate(kappa(p),p));
 286
 287 (kappa(x):=x, ev(%,nouns));
 288 V = %e^(p0^2/2-p^2/2)*V0;
 289
 290 /* Bug report ID: 1621 Wrong solution to ode2
 291  *
 292  * 'diff(y,t,2)-a*'diff(y,t)=-t with a=0
 293  *
 294  * Test case is a 2nd order non-homogeneous linear ode with constant
 295  * coefficients.  When solving homogeneous ode, cc2 asked if a=0, then
 296  * returned a solution involving a.
 297  */
 298 assume(equal(a,0));
 299 [equal(a,0)];
 300
 301 ode2('diff(y,t,2)-a*'diff(y,t)=0,y,t);
 302 y=%k2*t+%k1;
 303
 304 ode2('diff(y,t,2)-a*'diff(y,t)=-t,y,t);
 305 y=(-t^3/6)+%k2*t+%k1;
 306
 307 forget(equal(a,0));
 308 [equal(a,0)];
 309