| blob | ce3c0d22c9ec8a6b6a1d6715fa3ec1a863448796 |
1 /* Original version of this file copyright 1999 by Michael Wester,
2 * and retrieved from http://www.math.unm.edu/~wester/demos/ODEs/problems.macsyma
3 * circa 2006-10-23.
4 *
5 * Released under the terms of the GNU General Public License, version 2,
6 * per message dated 2007-06-03 from Michael Wester to Robert Dodier
7 * (contained in the file wester-gpl-permission-message.txt).
8 *
9 * See: "A Critique of the Mathematical Abilities of CA Systems"
10 * by Michael Wester, pp 25--60 in
11 * "Computer Algebra Systems: A Practical Guide", edited by Michael J. Wester
12 * and published by John Wiley and Sons, Chichester, United Kingdom, 1999.
13 */
14 /* ----------[ M a c s y m a ]---------- */
15 /* ---------- Initialization ---------- */
16 showtime: all$
17 prederror: false$
18 /* ---------- Ordinary Difference and Differential Equations ---------- */
19 /* Second order linear recurrence equation: r(n) = (n - 1)^2 + m n */
20 differenceq(r[n + 2] - 2 * r[n + 1] + r[n] = 2, r[n], [r[0] = 1, r[1] = m]);
21 /* => r(n) = 3^n - 2^n [Cohen, p. 67] */
22 errcatch(differenceq(r[n] = 5*r[n - 1] - 6*r[n - 2], r[n],
23 [r[0] = 0, r[1] = 1]))$
24 differenceq(r[n + 2] = 5*r[n + 1] - 6*r[n], r[n], [r[0] = 0, r[1] = 1]);
25 /* => r(n) = Fibonacci[n + 1] [Cohen, p. 83] */
26 errcatch(differenceq(r[n] = r[n - 1] + r[n - 2], r[n], [r[1] = 1, r[2] = 2]))$
27 expr: differenceq(r[n + 2] = r[n + 1] + r[n], r[n], [r[1] = 1, r[2] = 2]);
28 solve(2 = subst(n = 2, rhs(%)));
29 expr: factorsum(subst(%, expr));
30 subst(%phi = (1 + sqrt(5))/2, fibtophi(fib(n+1)));
31 ratsimp(rhs(expr) - %);
32 remvalue(expr)$
33 /* => [c^(n+1) [c^(n+1) - 2 c - 2] + (n+1) c^2 + 2 c - n] / [(c-1)^3 (c+1)]
34 [Joan Z. Yu and Robert Israel in sci.math.symbolic] */
35 eqn: r[n] = (1 + c - c^(n-1) - c^(n+1))/(1 - c^n)*r[n - 1]
36 - c*(1 - c^(n-2))/(1 - c^(n-1))*r[n - 2] + 1;
37 errcatch(differenceq(eqn, r[n], [r[1] = 1, r[2] = (2 + 2*c + c^2)/(1 + c)]))$
38 remvalue(eqn)$
39 /* Second order ODE with initial conditions---solve first using Laplace
40 transforms: f(t) = sin(2 t)/8 - t cos(2 t)/4 */
41 atvalue(f(t), t = 0, 0)$
42 atvalue(diff(f(t), t), t = 0, 0)$
43 printprops(all, atvalue)$
44 ode: diff(f(t), t, 2) + 4*f(t) = sin(2*t);
45 laplace(ode, t, s);
46 solve(%, 'laplace(f(t), t, s));
47 ilt(%[1], s, t);
48 /* Now, solve the ODE directly */
49 subst(f(t) = f, ode);
50 ode(%, f, t);
51 ode_ibc(%, t = 0, f = 0, t = 0, 'diff(f, t) = 0);
52 remvalue(ode)$
53 /* Separable equation => y(x)^2 = 2 log(x + 1) + (4 x + 3)/(x + 1)^2 + 2 A */
54 'diff(y, x) = x^2/(y*(1 + x)^3);
55 ode(%, y, x);
56 lhs(%) = map('factorsum, rhs(%));
57 /* Homogeneous equation. See Emilio O. Roxin, _Ordinary Differential
58 Equations_, Wadsworth Publishing Company, 1972, p. 11
59 => y(x)^2 = 2 x^2 log|A x| */
60 'diff(y, x) = y/x + x/y;
61 ode(%, y, x);
62 /* First order linear ODE: y(x) = [A - cos(x)]/x^3 */
63 x^2*'diff(y, x) + 3*x*y = sin(x)/x;
64 ode(%, y, x);
65 /* Exact equation => x + x^2 sin y(x) + y(x) = A [Roxin, p. 15] */
66 'diff(y, x) = -(1 + 2*x*sin(y))/(1 + x^2*cos(y));
67 ode(%, y, x);
68 /* Nonlinear ODE => y(x)^3/6 + A y(x) = x + B */
69 'diff(y, x, 2) + y*'diff(y, x)^3 = 0;
70 ode(%, y, x);
71 /* => y(x) = [3 x + sqrt(1 + 9 x^2)]^(1/3) - 1/[3 x + sqrt(1 + 9 x^2)]^(1/3)
72 [Pos96] */
73 errcatch(ode_ibc(%, x = 0, y = 0, x = 0, 'diff(y, x) = 2));
74 /* A simple parametric ODE: y(x, a) = A e^(a x) */
75 diff(y(x, a), x) = a*y(x, a);
76 ode(%, y(x, a), x);
77 /* ODE with boundary conditions. This problem has nontrivial solutions
78 y(x) = A sin([pi/2 + n pi] x) for n an arbitrary integer */
79 assume(not(equal(k, 0)))$
80 ode('diff(y, x, 2) + k^2*y = 0, y, x);
81 ode_ibc(%, x = 0, y = 0, x = 1, 'diff(y, x) = 0);
82 forget(not(equal(k, 0)))$
83 /* => y(x) = Z_v[sqrt(x)] where Z_v is an arbitrary Bessel function of order v
84 [Gradshteyn and Ryzhik 8.491(9)] */
85 eqn: 'diff(y, x, 2) + 1/x*'diff(y, x) + 1/(4*x)*(1 - v^2/x)*y = 0;
86 declare(v, noninteger)$
87 ode(eqn, y, x);
88 remove(v, noninteger)$
89 declare(v, integer)$
90 ode(eqn, y, x);
91 remove(v, integer)$
92 /* Delay (or mixed differential-difference) equation. See Daniel Zwillinger,
93 _Handbook of Differential Equations_, Second Edition, Academic Press, Inc.,
94 1992, p. 210 => y(t) = y0 sum((-a)^n (t - n + 1)^n/n!, n = 0..floor(t) + 1)
95 */
96 'diff(y(t), t) + a*y(t - 1) = 0;
97 ode(%, y(t), t);
98 /* Discontinuous ODE [Zwillinger, p. 221]
99 => y(t) = cosh t (0 <= t < T)
100 (sin T cosh T + cos T sinh T) sin t
101 + (cos T cosh T - sin T sinh T) cos t (T <= t) */
102 sgn(t):= if t < 0 then -1 else 1$
103 ode(diff(y(t), t, 2) + sgn(t - TT)*y(t) = 0, y(t), t);
104 /*ode_ibc(%, t = 0, y(t) = 1, t = 0, 'diff(y(t), t) = 0);*/
105 remfunction(sgn)$
106 assume(pnz > 0)$
107 ode(diff(y(t), t, 2) + sign(t - TT)*y(t) = 0, y(t), t);
108 errcatch(ode_ibc(%, t = 0, y(t) = 1, t = 0, 'diff(y(t), t) = 0));
109 forget(pnz > 0)$
110 /* Integro-differential equation. See A. E. Fitzgerald, David E. Higginbotham
111 and Arvin Grabel, _Basic Electrical Engineering_, Fourth Edition,
112 McGraw-Hill Book Company, 1975, p. 117.
113 => i(t) = 5/13 [-8 e^(-4 t) + e^(-t) (8 cos 2 t + sin 2 t)] */
114 eqn: diff(i(t), t) + 2*i(t) + 5*integrate(i(tau), tau, 0, t) = 10*%e^(-4*t);
115 ode(eqn, i(t), t);
116 atvalue(i(t), t = 0, 0)$
117 atvalue(diff(i(t), t), t = 0, 10)$
118 laplace(eqn, t, s);
119 solve(%, 'laplace(i(t), t, s));
120 ilt(%[1], s, t);
121 remvalue(eqn)$
122 /* System of two linear, constant coefficient ODEs:
123 x(t) = e^t [A cos(t) - B sin(t)], y(t) = e^t [A sin(t) + B cos(t)] */
124 system: [diff(x(t), t) = x(t) - y(t), diff(y(t), t) = x(t) + y(t)];
125 expand(odelinsys(system, [x(t), y(t)]));
126 /* Check the answer */
127 ratsimp(ev(system, %, diff));
128 map(lambda([q], is(equal(lhs(q), rhs(q)))), %);
129 /* Triangular system of two ODEs: x(t) = A e^t [sin(t) + 2],
130 y(t) = A e^t [5 - cos(t) + 2 sin(t)]/5 + B e^(-t).
131 See Nicolas Robidoux, ``Does Axiom Solve Systems of O.D.E.'s Like
132 Mathematica?'', LA-UR-93-2235, Los Alamos National Laboratory, Los Alamos,
133 New Mexico. */
134 system: [diff(x(t), t) = x(t) * (1 + cos(t)/(2 + sin(t))),
135 diff(y(t), t) = x(t) - y(t)];
136 odelinsys(system, [x(t), y(t)]);
137 /* Try solving this system one equation at a time */
138 subst(%c = %c1, factor(ode(system[1], x(t), t)));
139 map('factorsum, ode(subst(%, system[2]), y(t), t));
140 /* 3 x 3 linear system with constant coefficients:
141 (1) real distinct characteristic roots (= 2, 1, 3) [Roxin, p. 109]
142 => x(t) = A e^(2 t), y(t) = B e^t + C e^(3 t),
143 z(t) = -A e^(2 t) - C e^(3 t) */
144 system: [diff(x(t), t) = 2*x(t),
145 diff(y(t), t) = -2*x(t) + y(t) - 2*z(t),
146 diff(z(t), t) = x(t) + 3*z(t)];
147 odelinsys(system, [x(t), y(t), z(t)]);
148 /* (2) complex characteristic roots (= 0, -1 +- sqrt(2) i) [Roxin, p. 111]
149 => x(t) = A + e^(-t)/3 [-(B + sqrt(2) C) cos(sqrt(2) t) +
150 (sqrt(2) B - C) sin(sqrt(2) t)],
151 y(t) = e^(-t) [B cos(sqrt(2) t) + C sin(sqrt(2) t)],
152 z(t) = e^(-t) [(-B + sqrt(2) C) cos(sqrt(2) t)
153 -(sqrt(2) B + C) sin(sqrt(2) t)] */
154 system: [diff(x(t), t) = y(t), diff(y(t), t) = z(t),
155 diff(z(t), t) = -3*y(t) - 2*z(t)];
156 odelinsys(system, [x(t), y(t), z(t)]);
157 ratsimp(%);
158 /* (3) multiple characteristic roots (= 2, 2, 2) [Roxin, p. 113]
159 => x(t) = e^(2 t) [A + C (1 + t)], y(t) = B e^(2 t),
160 z(t) = e^(2 t) [A + C t] */
161 system: [diff(x(t), t) = 3*x(t) - z(t), diff(y(t), t) = 2*y(t),
162 diff(z(t), t) = x(t) + z(t)];
163 odelinsys(system, [x(t), y(t), z(t)]);
164 /* x(t) = x0 + [4 sin(w t)/w - 3 t] x0' [Rick Niles]
165 + 6 [w t - sin(w t)] y0 + 2/w [1 - cos(w t)] y0',
166 y(t) = -2/w [1 - cos(w t)] x0' + [4 - 3 cos(w t)] y0 + sin(w t)/w y0' */
167 system: [diff(x(t), t, 2) = 2*w*diff(y(t), t),
168 diff(y(t), t, 2) = -2*w*diff(x(t), t) + 3*w^2*y(t)];
169 odelinsys(system, [x(t), y(t)]);
170 remvalue(system)$