Mon Apr 6 20:12:04 MDT 1998 aquarius% mupad *----* MuPAD 1.4.0 -- Multi Processing Algebra Data Tool /| /| *----* | Copyright (c) 1997 - 98 by SciFace Software GmbH | *--|-* All rights reserved. |/ |/ *----* Licensed to: Michael Wester >> # ----------[ M u P A D ]---------- # >> # ---------- Initialization ---------- # >> TEXTWIDTH:= 80: >> read("../../Time.mupad"): >> # ---------- Ordinary Difference and Differential Equations ---------- # >> # Second order linear recurrence equationn: r(n) = (n - 1)^2 + m n # >> r(n + 2) - 2 * r(n + 1) + r(n) = 2; r(n) - 2 r(n + 1) + r(n + 2) = 2 Time: 90 msec Type: "_equal" >> solve(rec(%, r(n), {r(0) = 1, r(1) = m})); 2 {n + n (m - 2) + 1} Time: 10950 msec Type: DOM_SET >> # => r(n) = 3^n - 2^n [Cohen, p. 67] # >> solve(rec(r(n) = 5*r(n - 1) - 6*r(n - 2), r(n), {r(0) = 0, r(1) = 1})); n n {3 - 2 } Time: 790 msec Type: DOM_SET >> # => r(n) = Fibonacci[n + 1] [Cohen, p. 83] # >> solve(rec(r(n) = r(n - 1) + r(n - 2), r(n), {r(1) = 1, r(2) = 2})); { / 1/2 \n / 1/2 \n } { 1/2 | 5 | | 5 | } { (5 - 3) | 1/2 - ---- | 2 | ---- + 1/2 | } { \ 2 / \ 2 / } { -------------------------- - ----------------- } { 1/2 1/2 } { 5 - 5 5 - 5 } Time: 1430 msec Type: DOM_SET >> # => [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] # >> solve(rec(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, r(n), &> {r(1) = 1, r(2) = (2 + 2*c + c^2)/(1 + c)})); / / n - 2 | | c r(n - 2) (- c + 1 ) solve| rec| r(n) = - -------------------------- + | | n - 1 \ \ - c + 1 n - 1 n + 1 r(n - 1) (c - c - c + 1) ---------------------------------- + 1 , r(n), n - c + 1 { 2 } \ \ { 2 c + c + 2 } | | { r(1) = 1, r(2) = ------------ } | | { c + 1 } | | / / Time: 1690 msec Type: "solve" >> # Second order ODE with initial conditions---solve first using Laplace &> transforms: f(t) = sin(2 t)/8 - t cos(2 t)/4 # >> Ode:= f''(t) + 4*f(t) = sin(2*t); 4 f(t) + D(D(f))(t) = sin(2 t) Time: 390 msec Type: "_equal" >> transform::laplace(Ode, t, s); 2 2 4 laplace(f(t), t, s) - D(f)(0) - s f(0) + s laplace(f(t), t, s) = ------ 2 s + 4 Time: 2930 msec Type: "_equal" >> subs(op(linsolve([%], [transform::laplace(f(t), t, s)])), &> {f(0) = 0, D(f)(0) = 0}); 2 laplace(f(t), t, s) = -------------- 2 4 8 s + s + 16 Time: 770 msec Type: "_equal" >> transform::ilaplace(%, s, t); sin(2 t) t cos(2 t) f(t) = -------- - ---------- 8 4 Time: 1070 msec Type: "_equal" >> # Now, solve the ODE directly # >> solve(ode({Ode, f(0) = 0, D(f)(0) = 0}, f(t))); { sin(2 t) t cos(2 t) } { -------- - ---------- } { 8 4 } Time: 10610 msec Type: DOM_SET >> Ode:= NIL: >> # Separable equation => y(x)^2 = 2 log(x + 1) + (4 x + 3)/(x + 1)^2 + 2 A # >> y'(x) = x^2/(y(x)*(1 + x)^3); 2 x D(y)(x) = ------------- 3 y(x) (x + 1) Time: 610 msec Type: "_equal" >> solve(ode(%, y(x))); 2 2 2 {RootOf(2 C3 - 4 x + 4 C3 x + y - 2 ln(x + 1) + 2 C3 x + 2 x y - 2 2 2 4 x ln(x + 1) + x y - 2 x ln(x + 1) - 3, y)} Time: 7740 msec Type: DOM_SET >> # Homogeneous equation. See Emilio O. Roxin, _Ordinary Differential &> Equations_, Wadsworth Publishing Company, 1972, p. 11 &> => y(x)^2 = 2 x^2 log|A x| # >> y'(x) = y(x)/x + x/y(x); x y(x) D(y)(x) = ---- + ---- y(x) x Time: 320 msec Type: "_equal" >> solve(ode(%, y(x))); { / / 2 / x \ \ \ } { x RootOf| | - v2 + 2 ln| -- | | + (0, -4 I PI, 4 I PI, ...), v2 | } { \ \ \ C4 / / / } Time: 17170 msec Type: DOM_SET >> # First order linear ODE: y(x) = [A - cos(x)]/x^3 # >> x^2*y'(x) + 3*x*y(x) = sin(x)/x; 2 sin(x) 3 x y(x) + x D(y)(x) = ------ x Time: 360 msec Type: "_equal" >> solve(ode(%, y(x))); { C5 + cos(x) } { - ----------- } { 3 } { x } Time: 1230 msec Type: DOM_SET >> # Exact equation => x + x^2 sin y(x) + y(x) = A [Roxin, p. 15] # >> y'(x) = -(1 + 2*x*sin(y(x)))/(1 + x^2*cos(y(x))); 2 x sin(y(x)) + 1 D(y)(x) = - ----------------- 2 x cos(y(x)) + 1 Time: 360 msec Type: "_equal" >> solve(ode(%, y(x))); / / 2 x sin(y(x)) + 1 \ \ solve| ode| diff(y(x), x) = - -----------------, y(x) | | | | 2 | | \ \ x cos(y(x)) + 1 / / Time: 11070 msec Type: "solve" >> # Nonlinear ODE => y(x)^3/6 + A y(x) = x + B # >> eqn:= y''(x) + y(x)*y'(x)^3 = 0; 3 D(D(y))(x) + y(x) D(y)(x) = 0 Time: 400 msec Type: "_equal" >> solve(ode(%, y(x))); 3 {RootOf(6 x + 6 C8 + 6 y C6 - y , y), C7} Time: 1810 msec Type: DOM_SET >> # => y(x) = [3 x + sqrt(1 + 9 x^2)]^(1/3) - 1/[3 x + sqrt(1 + 9 x^2)]^(1/3) &> [Pos96] # >> solve(ode({eqn, y(0) = 0, D(y)(0) = 2}, y(x))); 3 solve(ode({y(0) = 0, D(y)(0) = 2, diff(y(x), x, x) + y(x) diff(y(x), x) = 0}, y(x))) Time: 2270 msec Type: "solve" >> # A simple parametric ODE: y(x, a) = A e^(a x) # >> diff(y(x, a), x) = a*y(x, a); diff(y(x, a), x) = a y(x, a) Time: 400 msec Type: "_equal" >> solve(ode(%, y(x, a))); {C12 exp(a x)} Time: 1310 msec Type: DOM_SET >> # ODE with boundary conditions. This problem has nontrivial solutions &> y(x) = A sin([pi/2 + n pi] x) for n an arbitrary integer # >> solve(ode({y''(x) + k^2*y(x) = 0, y(0) = 0, D(y)(1) = 0}, y(x))); {0} Time: 1500 msec Type: DOM_SET >> # => y(x) = Z_v[sqrt(x)] where Z_v is an arbitrary Bessel function of order v &> [Gradshteyn and Ryzhik 8.491(9)] # >> y''(x) + 1/x*y'(x) + 1/(4*x)*(1 - v^2/x)*y(x) = 0; / 2 \ | v | y(x) | 1 - -- | D(y)(x) \ x / D(D(y))(x) + ------- + --------------- = 0 x 4 x Time: 420 msec Type: "_equal" >> solve(ode(%, y(x))); / / / 2 \ \ \ | | | v | | | | | y(x) | - -- + 1 | | | | | diff(y(x), x) \ x / | | solve| ode| diff(y(x), x, x) + ------------- + ------------------ = 0, y(x) | | \ \ x 4 x / / Time: 26780 msec Type: "solve" >> # 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) # >> y'(t) + a*y(t - 1) = 0; D(y)(t) + a y(t - 1) = 0 Time: 450 msec Type: "_equal" >> solve(ode(%, y(t))); Error: not an ordinary differential equation in the given variable [ode::solve_eq_irred] >> # 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:= proc(t) begin if t < 0 then -1 else 1 end_if end_proc: >> solve(ode({y''(t) + sgn(t - TT)*y(t) = 0, y(0) = 1, D(y)(0) = 0}, y(t))); Error: Can't evaluate to boolean [_less]; in procedure 'sgn' >> sgn:= NIL: >> solve(ode({y''(t) + sign(t - TT)*y(t) = 0, y(0) = 1, D(y)(0) = 0}, y(t))); solve(ode({D(y)(0) = 0, y(0) = 1, diff(y(t), t, t) + y(t) sign(t - TT) = 0}, y(t))) Time: 20920 msec Type: "solve" >> # 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:= i'(t) + 2*i(t) + 5*int(i(tau), tau = 0..t) = 10*E^(-4*t); 10 2 i(t) + D(i)(t) + 5 int(i(tau), tau = 0..t) = --------- 4 t exp(1) Time: 4390 msec Type: "_equal" >> solve(ode(eqn, i(t))); { / 10 2 \ } { int| ------- - 5 exp(t) int(i(tau), tau = 0..t), t | } { | 2 | } { \ exp(t) / } { C16 exp(-2 t) + ----------------------------------------------------- } { 2 } { exp(t) } Time: 11930 msec Type: DOM_SET >> transform::laplace(eqn, t, s); 5 laplace(i(t), t, s) 2 laplace(i(t), t, s) + s laplace(i(t), t, s) - i(0) + --------------------- = s / 1 \ 10 laplace| ---------, t, s | | 4 t | \ exp(1) / Time: 1250 msec Type: "_equal" >> subs(op(linsolve([%], [transform::laplace(i(t), t, s)])), &> {i(0) = 0, D(i)(0) = 10}); Error: unable to convert i(0) + laplace(exp(1)^(t*(-4)), t, s)*10 [linalg::expr2Matrix] >> #transform::ilaplace(%, s, t);# >> eqn:= NIL: >> # 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:= {x'(t) = x(t) - y(t), y'(t) = x(t) + y(t)}; {D(x)(t) = x(t) - y(t), D(y)(t) = x(t) + y(t)} Time: 1560 msec Type: DOM_SET >> solve(ode(SysteM, {x(t), y(t)})); {{y(t) = C17 exp((1 - I) t) + C18 exp((1 + I) t), x(t) = I C18 exp((1 + I) t) - I C17 exp((1 - I) t)}} Time: 3130 msec Type: DOM_SET >> ans:= map(op(%), func(op(q, 1) = Factor(rectform(op(q, 2), indets(q))), q)); {y(t) = -exp(t) (I C17 sin(t) - C17 cos(t) - C18 cos(t) - I C18 sin(t)), x(t) = -exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t))} Time: 1350 msec Type: DOM_SET >> # Check the answer # >> subs(SysteM, op(ans)); {D(x)(t) = exp(t) (I C17 sin(t) - C17 cos(t) - C18 cos(t) - I C18 sin(t)) - exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t)), D(y)(t) = - exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t)) - exp(t) (I C17 sin(t) - C17 cos(t) - C18 cos(t) - I C18 sin(t))} Time: 540 msec Type: DOM_SET >> eval(subs({diff(x(t), t) = x(t) - y(t), diff(y(t), t) = x(t) + y(t)}, &> op(ans))); {- exp(t) (C17 cos(t) + C18 cos(t) - I C17 sin(t) + I C18 sin(t)) - exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t)) = exp(t) (I C17 sin(t) - C17 cos(t) - C18 cos(t) - I C18 sin(t)) - exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t)), - exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t)) - exp(t) (I C17 sin(t) - C17 cos(t) - C18 cos(t) - I C18 sin(t)) = - exp(t) (I C17 cos(t) - I C18 cos(t) + C17 sin(t) + C18 sin(t)) - exp(t) (I C17 sin(t) - C17 cos(t) - C18 cos(t) - I C18 sin(t))} Time: 580 msec Type: DOM_SET >> simplify(%, logic); {TRUE, FALSE} Time: 1680 msec Type: DOM_SET >> ans:= NIL: >> # 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:= [x'(t) = x(t) * (1 + cos(t)/(2 + sin(t))), &> y'(t) = x(t) - y(t)]; -- / cos(t) \ -- | D(x)(t) = x(t) | ---------- + 1 |, D(y)(t) = x(t) - y(t) | -- \ sin(t) + 2 / -- Time: 1040 msec Type: DOM_LIST >> solve(ode(SysteM, {x(t), y(t)})); Error: Recursive definition [See ?MAXDEPTH]; in procedure 'intlib::match_aux' >> # Try solving this system one equation at a time # >> x(t) = op(solve(ode(SysteM[1], x(t)))); 1/2 x(t) = C21 exp(t) (16 sin(t) - 2 cos(2 t) + 18) Time: 152520 msec Type: "_equal" >> subs(expand(%), cos(t)^2 = 1 - sin(t)^2); 2 1/2 x(t) = C21 exp(t) (16 sin(t) + 4 sin(t) + 16) Time: 580 msec Type: "_equal" >> map(%, func(Factor(q^2), q)); 2 2 2 2 x(t) = 4 C21 exp(t) (sin(t) + 2) Time: 840 msec Type: "_equal" >> x(t) = A*exp(t)*(sin(t) + 2); x(t) = A exp(t) (sin(t) + 2) Time: 570 msec Type: "_equal" >> y(t) = op(solve(ode(subs(SysteM[2], %), y(t)))); A cos(t) exp(t) 2 A sin(t) exp(t) y(t) = A exp(t) - --------------- + ----------------- + C22 exp(-t) 5 5 Time: 9440 msec Type: "_equal" >> # 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:= {x'(t) = 2*x(t), &> y'(t) = -2*x(t) + y(t) - 2*z(t), &> z'(t) = x(t) + 3*z(t)}; {D(x)(t) = 2 x(t), D(z)(t) = x(t) + 3 z(t), D(y)(t) = y(t) - 2 x(t) - 2 z(t)} Time: 590 msec Type: DOM_SET >> solve(ode(SysteM, {x(t), y(t), z(t)})); {{y(t) = C23 exp(t) - C25 exp(3 t), x(t) = -C24 exp(2 t), z(t) = C24 exp(2 t) + C25 exp(3 t)}} Time: 2730 msec Type: DOM_SET >> # (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:= {x'(t) = y(t), y'(t) = z(t), z'(t) = -3*y(t) - 2*z(t)}; {D(y)(t) = z(t), D(z)(t) = - 3 y(t) - 2 z(t), D(x)(t) = y(t)} Time: 610 msec Type: DOM_SET >> solve(ode(SysteM, {x(t), y(t), z(t)})); { { 1/2 { { 1/2 1/2 C27 exp(- t - I t 2 ) { { y(t) = 1/3 I C27 2 exp(- t - I t 2 ) - ----------------------- - { { 3 1/2 C28 exp(I t 2 - t) 1/2 1/2 --------------------- - 1/3 I C28 2 exp(I t 2 - t), 3 1/2 1/2 z(t) = C27 exp(- t - I t 2 ) + C28 exp(I t 2 - t), x(t) = C26 - 1/2 1/2 C27 exp(- t - I t 2 ) C28 exp(I t 2 - t) ----------------------- - --------------------- - 9 9 } } 1/2 1/2 1/2 1/2 } } 2/9 I C27 2 exp(- t - I t 2 ) + 2/9 I C28 2 exp(I t 2 - t) } } } } Time: 4620 msec Type: DOM_SET >> map(op(%), func(op(q, 1) = Factor(rectform(op(q, 2), indets(q))), q)); { { 1/2 1/2 1/2 { y(t) = (exp(-t) (I C27 sin(t 2 ) - C27 cos(t 2 ) - C28 cos(t 2 ) - { 1/2 1/2 1/2 1/2 1/2 I C28 sin(t 2 ) + I C27 2 cos(t 2 ) - I C28 2 cos(t 2 ) + 1/2 1/2 1/2 1/2 C27 2 sin(t 2 ) + C28 2 sin(t 2 ))) / 3, z(t) = - exp(-t) 1/2 1/2 1/2 1/2 (I C27 sin(t 2 ) - C27 cos(t 2 ) - C28 cos(t 2 ) - I C28 sin(t 2 )), 1/2 1/2 C27 exp(-t) cos(t 2 ) C28 exp(-t) cos(t 2 ) x(t) = C26 - ----------------------- - ----------------------- + 9 9 1/2 1/2 1/9 I C27 exp(-t) sin(t 2 ) - 1/9 I C28 exp(-t) sin(t 2 ) - 1/2 1/2 1/2 1/2 2/9 I C27 2 exp(-t) cos(t 2 ) + 2/9 I C28 2 exp(-t) cos(t 2 ) - 1/2 1/2 1/2 1/2 } 2 C27 2 exp(-t) sin(t 2 ) 2 C28 2 exp(-t) sin(t 2 ) } ------------------------------ - ------------------------------ } 9 9 } Time: 3400 msec Type: DOM_SET >> # (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:= {x'(t) = 3*x(t) - z(t), y'(t) = 2*y(t), z'(t) = x(t) + z(t)}; {D(y)(t) = 2 y(t), D(z)(t) = x(t) + z(t), D(x)(t) = 3 x(t) - z(t)} Time: 630 msec Type: DOM_SET >> solve(ode(SysteM, {x(t), y(t), z(t)})); {{z(t) = 0, y(t) = C29 exp(2 t), x(t) = 0}} Time: 1960 msec Type: DOM_SET >> # 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:= {x''(t) = 2*w*y'(t), &> y''(t) = -2*w*x'(t) + 3*w^2*y(t)}; 2 {D(D(x))(t) = 2 w D(y)(t), D(D(y))(t) = 3 w y(t) - 2 w D(x)(t)} Time: 630 msec Type: DOM_SET >> solve(ode(SysteM, {x(t), y(t)})); { { I C32 exp(-I t w) - I C33 exp(I t w) { { y(t) = ------------------------------------, { { w w C34 - 2 C32 exp(-I t w) - 2 C33 exp(I t w) } } x(t) = -------------------------------------------- } } w } } Time: 5040 msec Type: DOM_SET >> map(op(%), func(op(q, 1) = Factor(rectform(op(q, 2), indets(q))), q)); { { x(t) = { w C34 - 2 C32 cos(t w) - 2 C33 cos(t w) + 2 I C32 sin(t w) - 2 I C33 sin(t w) ----------------------------------------------------------------------------- w I C32 cos(t w) - I C33 cos(t w) + C32 sin(t w) + C33 sin(t w) } , y(t) = ------------------------------------------------------------- } w } Time: 1720 msec Type: DOM_SET >> SysteM:= NIL: >> # ---------- Quit ---------- # >> quit real 346.66 user 342.55 sys 1.64