Solving Ordinary Differential Equations with Maple...
Getting started...
Type maple to begin.
oscar.Princeton.EDU% maple
|\^/| Maple V Release 5 (WMI Campus Wide License)
._|\| |/|_. Copyright (c) 1981-1997 by Waterloo Maple Inc. All rights
\ MAPLE / reserved. Maple and Maple V are registered trademarks of
<____ ____> Waterloo Maple Inc.
| Type ? for help.
First Order Linear ODEs
# y'(x) + (1/x) * y(x) = 1
> sol1 := dsolve(diff(y(x), x) + y(x) / x = 1, y(x));
_C1
sol1 := y(x) = 1/2 x + ---
x
#This is a general solution
# Let's apply an initial condition y(1) = -1 and find the constant _C1
> dsolve({diff(y(x), x) + y(x) / x =1 , y(1) = -1} , y(x));
y(x) = 1/2 x - 3/2 1/x
# Thus _C1 = -3/2
# Another example
# y'(x) = 8 * x^3 * y^2
> dsolve(diff(y(x), x) = 8 * x^3 * y(x)^2, y(x));
1
y(x) = - ----------
4
2 x - _C1
Second Order ODEs
# y''(x) + y(x) = 0
# Notice that $2 indicates a second derivative
> eq1 := diff(y(x), x$2) + y(x) = 0;
/ 2 \
|d |
eq1 := |--- y(x)| + y(x) = 0
| 2 |
\dx /
#Let's find a general solution
> sol1 := dsolve(eq1, y(x));
sol1 := y(x) = _C1 sin(x) + _C2 cos(x)
# Note two constants: _C1 and _C2
# Let's try another example...
# y''(x) + y(x) = x
> eq2 := diff(y(x), x$2) + y(x) = x;
/ 2 \
|d |
eq2 := |--- y(x)| + y(x) = x
| 2 |
\dx /
> sol2 := dsolve(eq2, y(x));
sol2 := y(x) = x + _C1 sin(x) + _C2 cos(x)
>
Using Initial Conditions...
# Consider the following equation with initial conditions:
# y'' + y = sin(t)
# y(0) = 0 and y'(0) = 1
> eq5 := dsolve({diff(y(t), t$2) + y(t) = sin(t), y(0) = 0, D(y)(0) = 1}, y(t));
3
eq5 := y(t) = 1/2 sin(t) + (1/2 cos(t) sin(t) - 1/2 t) cos(t) + sin(t)
# Notice that there are no arbitrary constants in this solution
# Function rhs() is used to obtain the right hand side of eq5 in the example below.
# Type this to plot the graph
> plot(rhs(eq5), t=0..4*Pi);
|
+ AAAA
| AA A
+ A A
4+ A A
| A A
+ A A
+ A A
| AA A
+ A A
+ A A
| A AA
2+ AAAAAA A A
+ AAAA AAA A A
| AAA AA A A
+ AA A A A
+ AA AA AA A
| AA A A A
+ AAA A A A
|AA 2 4A 6 A 8 10 A 12
*--+--+--+--+--+--+--+--+--+--+*-+--+--+--+--+--+--+--+*-+--+--+--+--+--+--+--+-*--+--+--+--+--
0+ A A AA
| A A A
+ A A A
+ AA A A
| A A A
+ A A A
-2+ A AA AA
| A A A
+ AA A A
+ AA A AA
| AA AA A
+ AAAAAA A
+ A
| A
-4+ A
+ A
| A
+ A
| A
+ A
+ A
| A
-6+ AA
+ AAAA
|
Solving Systems of ODEs
# Consider the system of ODEs
x' = 4x + 2y
y' = 3x + 3y
> ex1 := {diff(x1(t), t) = 4 * x1(t) + 2 * x2(t),
> diff(x2(t), t) = 3 * x1(t) + 3 * x2(t)};
d d
ex1 := {-- x2(t) = 3 x1(t) + 3 x2(t), -- x1(t) = 4 x1(t) + 2 x2(t)}
dt dt
>
> dsolve(ex1, {x1(t), x2(t)} );
{x1(t) = 2/5 _C1 exp(t) + 3/5 _C1 exp(6 t) + 2/5 _C2 exp(6 t) - 2/5 _C2 exp(t),
x2(t) = 3/5 _C1 exp(6 t) - 3/5 _C1 exp(t) + 3/5 _C2 exp(t) + 2/5 _C2 exp(6 t) }
Numerical Solutions of ODEs
# Some problems do not have solutions in closed-form.
# Consider the motion of a simple pendulum with unit natural frequency:
# y'' + sin(y) = 0
> pend := diff(y(x), x$2) + sin(y(x)) = 0;
/ 2 \
|d |
pend := |--- y(x)| + sin(y(x)) = 0
| 2 |
\dx /
>
# Let's try to solve it...
> dsolve(pend, y(x));
y(x) y(x)
/ /
| 1 | 1
| -------------------- d_a - x - _C2 = 0, - | -------------------- d_a - x - _C2 = 0
| 1/2 | 1/2
/ (2 cos(_a) + _C1) / (2 cos(_a) + _C1)
# Maple found two implicit solutions, expressing 'x' as a function of 'y' in terms of
# elliptic integrals.
# Let's find the numerical solution to the pendulum equations.
# Suppose that y(0) = 0 and y'(0) = 1.
> sol := dsolve( {pend, y(0) = 0, D(y)(0) = 1}, y(x), type=numeric);
sol := proc(rkf45_x) ... end
# Note that the solution is returned as a procedure rkf45_x, displayed in abbreviated form.
# We will use function odeplot() to get a sketch of the result.
# This function is given here with tree arguments:
# - the name of the procedure returned by dsolve(),
# - a list with the names of the independent and dependent variables, and
# - a range for the independent variable.
> plots[odeplot](sol, [x, y(x)], 0..3*Pi);
|
| AAAAA AAAA
1+ AA A AA AA
| A A A A
| AA A AA A
+ A A A AA
| A A A A
| A A AA A
+ A A A A
| AA A AA A
| A AA A AA
+ A A A A
| A A A A
| A A A
+ A A A
| A A A
0.5+ A A A
| A A A
| A A A
+ A A A
| A A A
| A A A
+ A A A
| A A A
| A A A
+ A AA A
| A A A
|AA A A
+A A A
|A A A
|A A A
--*---+---+---+--+---+---+---+---+-*+---+---+---+---+--+---+---+--*+---+--+---+---+---+---+--+-
0| 2 A 4 6 A 8
+ A A
| A A
| A A
+ A A
| A A
| A A
+ A A
| A A
| A A
+ A A
| A A
| A A
-0.5+ A A
| A A
| A A
+ A A
| A A
+ A A
| A AA
| A A
+ A A
| A A
| A A
+ A A
| AA A
| AA AA
-1+ AA AA
| AAAA
|
>
>
Numerical Solutions of Systems of ODEs
# y' = x
# x' = -y - 0.1 * x * (y^2-1)
# suppose y(0) = 1/2 and x(0) = 0
> eq9 := {diff(y(t), t) = x(t),
> diff(x(t), t) = -y(t) - 0.1 * (y(t)^2 - 1) * x(t),
> y(0) = 1/2, x(0) = 0};
d d 2
eq9 := {-- y(t) = x(t), -- x(t) = -y(t) - .1 (y(t) - 1) x(t), y(0) = 1/2, x(0) = 0}
dt dt
>
> sol5 := dsolve(eq9, {x(t), y(t)}, type = numeric);
sol5 := proc(rkf45_x) ... end
>
>
> plots[odeplot](sol5, [[t, x(t)], [t, y(t)]], 0..20);
|
+ BB
| AAA B B
1+ A A B B
| A A B B
+ A A B B
| BBBB A AB B
+ AAA B B A * B
| AA A BB B AA * B
+ A A B BB A *A B
| BBB A AB B A BA B
+ AAA B B A BA B A B A B
| A AA B BB A BA B A B A B
0.5*B A AB B A BBA B A B A B
| BB A *A B A B A B A B A B
+ B AA BA B A B A B A B A
| B A B A BB A B A B A B A
+ BB A B A B A B A B A B AA
| B A BB A B A B A B A B A
+ B AA B A B A B A B AA B A
| B A B A B A B A B A B A
+ B A B AA B AA B A B A B A
| B A B A B A B A B A B A
--**---+-*-+---+*---+--*+----+*--+---*+---+-*-+----*---+--*-+---+*---+--*+---+*---+---*----+--
0|A B AA B 5 A B A10 B A B 15 A B A 20
+A B A B A B A B A B A B A
| A B A B A B A B A B A B A
+ A B A B A B AA B A B A B A
| A B A B A B A B A B A B A
+ A BA B A B A B A B A B A
| A BA B A B A B A B A B A
+ A AB B A BA B A B A B A
-0.5+ AA A BB BB AA BA B A B A BB A
| AAA B BB A *B B A B A B A
+ BBB A AB B A BA B A
| A AA B BB AA * B A
+ AAA B B A * B A
| A BBBB A AB B A
+ B A AB BB A
| A A BB B A
+ AAA B BB A
| A B B A
-1+ BBB A
| A
|
>
Displaying Graphs in Separate Windows
> plotsetup(x11,xterm);
> plot(x^2, x=-2..2);
These examples were taken from [Robertson 1996].
Yefim Shuf |
CS Department |
Princeton University
Last modified: Wed Oct 21 13:20:24 EDT 1998