/* ----------[ M a c s y m a ]---------- */ /* ---------- Initialization ---------- */ showtime: all$ prederror: false$ /* ---------- Ordinary Difference and Differential Equations ---------- */ /* Second order linear recurrence equation: r(n) = (n - 1)^2 + m n */ differenceq(r[n + 2] - 2 * r[n + 1] + r[n] = 2, r[n], [r[0] = 1, r[1] = m]); /* => r(n) = 3^n - 2^n [Cohen, p. 67] */ errcatch(differenceq(r[n] = 5*r[n - 1] - 6*r[n - 2], r[n], [r[0] = 0, r[1] = 1]))$ differenceq(r[n + 2] = 5*r[n + 1] - 6*r[n], r[n], [r[0] = 0, r[1] = 1]); /* => r(n) = Fibonacci[n + 1] [Cohen, p. 83] */ errcatch(differenceq(r[n] = r[n - 1] + r[n - 2], r[n], [r[1] = 1, r[2] = 2]))$ expr: differenceq(r[n + 2] = r[n + 1] + r[n], r[n], [r[1] = 1, r[2] = 2]); solve(2 = subst(n = 2, rhs(%))); expr: factorsum(subst(%, expr)); subst(%phi = (1 + sqrt(5))/2, fibtophi(fib(n+1))); ratsimp(rhs(expr) - %); remvalue(expr)$ /* => [c^(n+1) [c^(n+1) - 2 c - 2] + (n+1) c^2 + 2 c - n] / [(c-1)^3 (c+1)] [Joan Z. Yu and Robert Israel in sci.math.symbolic] */ eqn: r[n] = (1 + c - c^(n-1) - c^(n+1))/(1 - c^n)*r[n - 1] - c*(1 - c^(n-2))/(1 - c^(n-1))*r[n - 2] + 1; errcatch(differenceq(eqn, r[n], [r[1] = 1, r[2] = (2 + 2*c + c^2)/(1 + c)]))$ remvalue(eqn)$ /* Second order ODE with initial conditions---solve first using Laplace transforms: f(t) = sin(2 t)/8 - t cos(2 t)/4 */ atvalue(f(t), t = 0, 0)$ atvalue(diff(f(t), t), t = 0, 0)$ printprops(all, atvalue)$ ode: diff(f(t), t, 2) + 4*f(t) = sin(2*t); laplace(ode, t, s); solve(%, 'laplace(f(t), t, s)); ilt(%[1], s, t); /* Now, solve the ODE directly */ subst(f(t) = f, ode); ode(%, f, t); ode_ibc(%, t = 0, f = 0, t = 0, 'diff(f, t) = 0); remvalue(ode)$ /* Separable equation => y(x)^2 = 2 log(x + 1) + (4 x + 3)/(x + 1)^2 + 2 A */ 'diff(y, x) = x^2/(y*(1 + x)^3); ode(%, y, x); lhs(%) = map('factorsum, rhs(%)); /* Homogeneous equation. See Emilio O. Roxin, _Ordinary Differential Equations_, Wadsworth Publishing Company, 1972, p. 11 => y(x)^2 = 2 x^2 log|A x| */ 'diff(y, x) = y/x + x/y; ode(%, y, x); /* First order linear ODE: y(x) = [A - cos(x)]/x^3 */ x^2*'diff(y, x) + 3*x*y = sin(x)/x; ode(%, y, x); /* Exact equation => x + x^2 sin y(x) + y(x) = A [Roxin, p. 15] */ 'diff(y, x) = -(1 + 2*x*sin(y))/(1 + x^2*cos(y)); ode(%, y, x); /* Nonlinear ODE => y(x)^3/6 + A y(x) = x + B */ 'diff(y, x, 2) + y*'diff(y, x)^3 = 0; ode(%, y, x); /* => y(x) = [3 x + sqrt(1 + 9 x^2)]^(1/3) - 1/[3 x + sqrt(1 + 9 x^2)]^(1/3) [Pos96] */ errcatch(ode_ibc(%, x = 0, y = 0, x = 0, 'diff(y, x) = 2)); /* A simple parametric ODE: y(x, a) = A e^(a x) */ diff(y(x, a), x) = a*y(x, a); ode(%, y(x, a), x); /* ODE with boundary conditions. This problem has nontrivial solutions y(x) = A sin([pi/2 + n pi] x) for n an arbitrary integer */ assume(not(equal(k, 0)))$ ode('diff(y, x, 2) + k^2*y = 0, y, x); ode_ibc(%, x = 0, y = 0, x = 1, 'diff(y, x) = 0); forget(not(equal(k, 0)))$ /* => y(x) = Z_v[sqrt(x)] where Z_v is an arbitrary Bessel function of order v [Gradshteyn and Ryzhik 8.491(9)] */ eqn: 'diff(y, x, 2) + 1/x*'diff(y, x) + 1/(4*x)*(1 - v^2/x)*y = 0; declare(v, noninteger)$ ode(eqn, y, x); remove(v, noninteger)$ declare(v, integer)$ ode(eqn, y, x); remove(v, integer)$ /* Delay (or mixed differential-difference) equation. See Daniel Zwillinger, _Handbook of Differential Equations_, Second Edition, Academic Press, Inc., 1992, p. 210 => y(t) = y0 sum((-a)^n (t - n + 1)^n/n!, n = 0..floor(t) + 1) */ 'diff(y(t), t) + a*y(t - 1) = 0; ode(%, y(t), t); /* Discontinuous ODE [Zwillinger, p. 221] => y(t) = cosh t (0 <= t < T) (sin T cosh T + cos T sinh T) sin t + (cos T cosh T - sin T sinh T) cos t (T <= t) */ sgn(t):= if t < 0 then -1 else 1$ ode(diff(y(t), t, 2) + sgn(t - TT)*y(t) = 0, y(t), t); /*ode_ibc(%, t = 0, y(t) = 1, t = 0, 'diff(y(t), t) = 0);*/ remfunction(sgn)$ assume(pnz > 0)$ ode(diff(y(t), t, 2) + sign(t - TT)*y(t) = 0, y(t), t); errcatch(ode_ibc(%, t = 0, y(t) = 1, t = 0, 'diff(y(t), t) = 0)); forget(pnz > 0)$ /* Integro-differential equation. See A. E. Fitzgerald, David E. Higginbotham and Arvin Grabel, _Basic Electrical Engineering_, Fourth Edition, McGraw-Hill Book Company, 1975, p. 117. => i(t) = 5/13 [-8 e^(-4 t) + e^(-t) (8 cos 2 t + sin 2 t)] */ eqn: diff(i(t), t) + 2*i(t) + 5*integrate(i(tau), tau, 0, t) = 10*%e^(-4*t); ode(eqn, i(t), t); atvalue(i(t), t = 0, 0)$ atvalue(diff(i(t), t), t = 0, 10)$ laplace(eqn, t, s); solve(%, 'laplace(i(t), t, s)); ilt(%[1], s, t); remvalue(eqn)$ /* System of two linear, constant coefficient ODEs: x(t) = e^t [A cos(t) - B sin(t)], y(t) = e^t [A sin(t) + B cos(t)] */ system: [diff(x(t), t) = x(t) - y(t), diff(y(t), t) = x(t) + y(t)]; expand(odelinsys(system, [x(t), y(t)])); /* Check the answer */ ratsimp(ev(system, %, diff)); map(lambda([q], is(equal(lhs(q), rhs(q)))), %); /* Triangular system of two ODEs: x(t) = A e^t [sin(t) + 2], y(t) = A e^t [5 - cos(t) + 2 sin(t)]/5 + B e^(-t). See Nicolas Robidoux, ``Does Axiom Solve Systems of O.D.E.'s Like Mathematica?'', LA-UR-93-2235, Los Alamos National Laboratory, Los Alamos, New Mexico. */ system: [diff(x(t), t) = x(t) * (1 + cos(t)/(2 + sin(t))), diff(y(t), t) = x(t) - y(t)]; odelinsys(system, [x(t), y(t)]); /* Try solving this system one equation at a time */ subst(%c = %c1, factor(ode(system[1], x(t), t))); map('factorsum, ode(subst(%, system[2]), y(t), t)); /* 3 x 3 linear system with constant coefficients: (1) real distinct characteristic roots (= 2, 1, 3) [Roxin, p. 109] => x(t) = A e^(2 t), y(t) = B e^t + C e^(3 t), z(t) = -A e^(2 t) - C e^(3 t) */ system: [diff(x(t), t) = 2*x(t), diff(y(t), t) = -2*x(t) + y(t) - 2*z(t), diff(z(t), t) = x(t) + 3*z(t)]; odelinsys(system, [x(t), y(t), z(t)]); /* (2) complex characteristic roots (= 0, -1 +- sqrt(2) i) [Roxin, p. 111] => x(t) = A + e^(-t)/3 [-(B + sqrt(2) C) cos(sqrt(2) t) + (sqrt(2) B - C) sin(sqrt(2) t)], y(t) = e^(-t) [B cos(sqrt(2) t) + C sin(sqrt(2) t)], z(t) = e^(-t) [(-B + sqrt(2) C) cos(sqrt(2) t) -(sqrt(2) B + C) sin(sqrt(2) t)] */ system: [diff(x(t), t) = y(t), diff(y(t), t) = z(t), diff(z(t), t) = -3*y(t) - 2*z(t)]; odelinsys(system, [x(t), y(t), z(t)]); ratsimp(%); /* (3) multiple characteristic roots (= 2, 2, 2) [Roxin, p. 113] => x(t) = e^(2 t) [A + C (1 + t)], y(t) = B e^(2 t), z(t) = e^(2 t) [A + C t] */ system: [diff(x(t), t) = 3*x(t) - z(t), diff(y(t), t) = 2*y(t), diff(z(t), t) = x(t) + z(t)]; odelinsys(system, [x(t), y(t), z(t)]); /* x(t) = x0 + [4 sin(w t)/w - 3 t] x0' [Rick Niles] + 6 [w t - sin(w t)] y0 + 2/w [1 - cos(w t)] y0', y(t) = -2/w [1 - cos(w t)] x0' + [4 - 3 cos(w t)] y0 + sin(w t)/w y0' */ system: [diff(x(t), t, 2) = 2*w*diff(y(t), t), diff(y(t), t, 2) = -2*w*diff(x(t), t) + 3*w^2*y(t)]; odelinsys(system, [x(t), y(t)]); remvalue(system)$ /* ---------- Quit ---------- */ quit();