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