Edwards-Penney Chapter 4, sections 4.1, 4.2, 4.3, 4.4 Topics,
Definitions, Theorems, Examples
4.1 First order systems and applications
====
EXAMPLE. Forced system of undamped oscillators
m1 x'' = - k1 x + k2 (y-x),
m2 y'' = -k2(y-x) + f(t)
INSTANCE m1=2, m2=1, k1=4, k2=2, f(t)=40sin(3t)
2 x'' = - 6x + 2y,
y'' = 2x - 2y + 40 sin(3t)
EXAMPLE. Recirculating brine tanks
20 x' = -6x + y,
20 y' = 6x - 3y
x(t)=pounds of salt in tank 1 (100 gal)
y(t)=pounds of salt in tank 2 (200 gal)
x(0), y(0) = initial salt amounts in each tank
t=minutes
20=inflow rate=outflow rate
0=inflow salt concentration
FIRST ORDER SYSTEMS
EXAMPLE. Convert x''' + 3x'' + 2x' - 5x = sin(2t) into a first order
system using the position-velocity substitution x1=x, x2=x',
x3=x''.
ANSWER: x1' = x2,
x2' = x3,
x3' = 5x1 - 2x2 - 3x3 + sin(2t)
EXAMPLE. Transform to a first order system x'' = x^3 + (x')^3
ANSWER: x1=x, x2=x'', x1'=x2, x2'=(x1)^3+(x2)^3
EXAMPLE. Transform to a first order system
2x'' = -6x + 2y,
y'' = 2x - 2y + 40 sin(3t)
ANSWER: x1=x,x2=x',x3=y,x4=y' ==>
x1' = x2,
x2' = -3x1 + x3, [a division by 2 needed]
x3' = x4,
x4' = 2x1 - 2x3 + 40 sin(3t)
DYNAMICAL SYSTEMS
The typical second order equation from circuit theory and mechanical
oscillations is x'' + px' + qx = 0. It is realized as a 2-dimensional
system
x' = y,
y' = -qx - py
EXAMPLE. Solve x'=-2y, y'=x/2.
ANSWER: x(t)=A cos(t) + B sin(t), y(t) = (-A/2) cos(t) +
(B/2) sin(t).
METHOD: Differentiate the first equation, substitute the
second to get x''+x=0. Then x(t) is a linear combination of
cos(t), sin(t). From the first equation, y=(-1/2)x', which
implies the equation for y(t) by differentiation.
This method is called the CAYLEY-HAMILTON SHORTCUT.
THEOREM. [Cayley-Hamilton-Ziebur]
A linear dynamical system u'=Au with 2x2 constant matrix A can be
solved by finding the roots of the characteristic equation
|A-rI|=0. Then x(t) and y(t) are linear combinations of the
solution atoms from Euler's Theorem corresponding to the roots.
EXAMPLE. Solve x'=-2y, y'=x/2, x(0)=2, y(0)=0.
ANSWER: x(t)=2 cos(t), y(t)=sin(t)
METHOD: Find A,B in the general solution from the two
initial conditions at t=0.
EXAMPLE. Solve x'=y, y'=2x+y
ANSWER: x(t)=A exp(-t) + B exp(2t),
y(t)=-a exp(-t)+2B exp(2t)
METHOD: Cayley-Hamilton shortcut; A:=Matrix([[0,1],[2,1]]).
Char equation r^2-r-2=0, roots r=-1,2, solution atoms are
exp(-t), exp(2t). Then x(t)=linear combination. To find y,
use the first equation, y=x'.
EXAMPLE. Solve x'=-y, y'=1.01 x - 0.2 y, x(0)=0, y(0)=1.
ANSWER: x(t)=exp(-t/10) sin(t),
y(t)=(1/10)exp(-t/10)(sin(t)+10 cos(t))
METHOD: Cayley-Hamilton-Ziebur shortcut;
A:=Matrix([[0,-1],[1.01,-0.2]]) with char equation
r^2+0.2r+1.01=0 having complex roots 0.1+i, 0.1-i and
corresponding Euler solution atoms exp(-t/10)cos(t),
exp(-t/10)sin(t). Then x(t) is a linear combination of the
solution atoms and y=-x' from the first equation. To find
A,B in the general solution, use the two initial conditions
at t=0.
THEOREM 1. [Picard's Existence-Uniqueness for linear systems]
Consider the system of differential equations u'=P(t)u+F(t) where
u=, P(t) is nxn with continuous coefficients,
F(t)= has continuous entries. Let b= be a
constant vector and t=a a point in the interval J of continuity. Then
the problem
u' = P(t)u + F(t),
u(a)=b,
has a unique solution u(t) defined on J.
4.2 The Method of Elimination [Cayley-Hamilton-Ziebur Method]
====
EXAMPLE 1. Solve x' = 4x - 3y, y' = 6x - 7y, x(0)=2, y(0)=-1.
SOLUTION: Elimination is equivalent to using Cayley-Hamilton-Ziebur.
This example shows why. The characteristic equation is
identical to the operational determinant, with roots
r=2,-5 and Euler solution atoms e^(2t), e^(-5t)$. Then
we have the conclusion of C-H_Z, that x(t) and y(t) are
linear combinations of the Euler atoms:
x(t) = a1 e^(2t) + a2 e^(-5t),
y(t) = b1 e^(2t) + b2 e^(-5t).
Exactly as in C-H-Z, these equations are substituted
into the first differential equation, then coefficients
of atoms are matched to get equations for a1,a2,b1,b2.
ANSWER: x(t) = (3/2) b1 e^(2t) + (1/3) b2 e^(-5t),
y(t) = b1 e^(2t) + b2 e^(-5t).
EXAMPLE 2. Solve x'' + 3x - y = 0, -2x + y'' + 2y = 0.
SOLUTION: The system can be represented as a 4x4 linear system
of the form Au=0. The operational determinant is
identical to |A - rI|, which is the characteristic
equation of matrix A.
According to Cayley-Hamilton-Ziebur, the components x,y
are linear combinations of the Euler atoms found from
the four roots i, -i, 2i, -2i of |A-rI|=0.
The shortcut applies again, using the equations
x(t) = a1 cos(t) + a2 sin(t) + a3 cos(2t) + a4 sin(2t),
y(t) = b1 cos(t) + b2 sin(t) + b3 cos(2t) + b4 sin(2t).
The plan is to substitute these equations into the first
differential equation, then match coefficients of atoms
to obtain formulas for b1 to b4 in terms of a1 to a4.
ANSWER: x(t) = a1 cos(t) + a2 sin(t) + a3 cos(2t) + a4 sin(2t),
y(t) = 2a1 cos(t) + 2a2 sin(t) - b1 cos(2t) - b2 sin(2t).
DEFINITION: The NATURAL FREQUENCIES of oscillation are the trig
frequencies 1, 2 in the example.
See Example 1, section 4.1.
4.3 Numerical Methods for Systems
====
DE: du/dt = F(t,u), u=vector, t=scalar
EULER METHOD.
The vector formulas for iteration are
t[n+1] = t[n] + h
u[n+1] = u[n] + h F(t[n],u[n])
Given are u[0], t[0], h. The symbol F is the vector function
defining the differential equation du/dt = F(t,u).
IMPROVED EULER METHOD. HEUN METHOD.
The vector formula for iteration is
t[n+1] = t[n] + h Time variable
y[n+1] = u[n] + h * F(t[n],u[n]) Predictor
u[n+1] = u[n] + h/2 * F(t[n],u[n])
+ h/2 * F(t[n+1],y[n+1]) Corrector
Given are u[0], t[0], h. The symbol F is the vector function
defining the differential equation du/dt = F(t,u). Intermediate
symbol y[n] is used only to pass the predictor value to the
corrector.
RK4 METHOD.
The vector formula for iteration has five lines:
t[n+1] = t[n] + h
k1 = F(t[n],u[n]) Predictor 1
k2 = F(t[n]+h/2,u[n] + k1*h/2) Predictor 2
k3 = F(t[n]+h/2,u[n] + k2*h/2) Predictor 3
k4 = F(t[n+1],u[n] + k3*h) Predictor 4
u[n+1] = u[n] + h(k1 + 2*k2 + 2*k3 + k4)/6 Corrector
Given are u[0], t[0], h. The symbol F is the vector function
defining the differential equation du/dt = F(t,u). Intermediate
predictors k1, k2, k3, k4 are only used for the corrector, the
fifth formula, their values forgotten after one use.
EXAMPLE 1. Solve numerically the system
du/dt = F(t,u), u[0] = (3,6)
where u=(x(t),y(t)), F(t,u) = (3x-2y,5x-4y).
Compare the calculations to the exact solution
x(t)=2 e^(-2t) + e^t, y(t) = 5e^(-2t) + e^t
using the Euler Method to find values for u(0.1),
u(0.2) and u(0.3).
ANSWERS: u[0]=(3,6), t[0]=0, h=0.1. A hand calculator
is useful. Then u(0.1) approx u[1], u(0.2) approx u[2]
and u(0.3) approx u[3].
u[1] = u[0] + h*F(0,u[0]) = (2.7,5.1)
u[2] = u[1] + h*F(h,u[1]) = (2.4,4.2)
u[3] = u[2] + h*F(2h,u[2]) = (2.58,4.62)
EXAMPLE 2. Solve numerically the system
du/dt = F(t,u), u[0] = (0,1)
where u=(x(t),y(t)), F(t,u) = (y,-x).
Compare the calculations to the exact solution
x(t)=sin(t), y(t) = cos(t).
Use the vector RK4 Method to find values for u(t)
from t=0 to t=5 in steps h=0.5.
ANSWER: See table 4.3.1 in the textbook.
EXAMPLE 3. Solve numerically the Lunar Lander system
du/dt = F(t,u), u[0] = (1781870,-450)
where u=(x(t),y(t)), F(t,u) = (y,4 - 4.9044*10^12 /x^2).
Use the RK4 Method with h=0.1.
COMPUTER CODE: See the section on comets and
spacecraft.
Maple code for an answer check:
des:=diff(x(t),t)=y(t),diff(y(t),t)=4 - 4.9044*10^12 /x(t)^2;
ic:=x(0)=1781870,y(0)=-450;
p:=dsolve({des,ic},{x(t),y(t)},numeric); p(100); p(200);
EXAMPLE 4. Baseball problem. The equations are
du/dt = F(t,u)
u(0)=u[0]=(0,0,80*sqrt(3),80)
h=0.1, 0.05 and 0.025 for multiple computations
F(t,u) = (f1,f2,f3,f4)
f1 = u3
f2 = u4
f3 = -c*u3*sqrt(u3^2+u4^2)
f4 = -g - c*u4*sqrt(u3^2+u4^2)
c = 0 No Air Resistance
c = 0.0025 Nonlinear Newton drag
g = 32
ANSWERS: For c=0 the exact solution is x(t)=80*sqrt(3)*t,
y(t)=80t-16*t^2. Here's a quick answer check:
des0:=diff(x(t),t,t)=0,diff(y(t),t,t)=-32;
ic:=x(0)=0,y(0)=0,D(x)(0)=80*sqrt(3),D(y)(0)=80;
dsolve({des0,ic},{x(t),y(t)});
For c=0.0025, no exact solution was attempted. Instead, we
code for maple:
c=0.0025;
des1:=diff(x(t),t,t)=-c*sqrt((D(x)(t))^2+(D(y)(t))^2)*D(x)(t),
diff(y(t),t,t)=-32-c*sqrt((D(x)(t))^2+(D(y)(t))^2)*D(y)(t);
ic:=x(0)=0,y(0)=0,D(x)(0)=80*sqrt(3),D(y)(0)=80;
p:=dsolve({des1,ic},{x(t),y(t)},numeric);
An answer check for Table 4.3.5 columns 1,2 can be made by
printing p(t), for t=0 to t=4 step 0.5, e.g., p(4.0) is
[t = 4.,
x(t) = 339.671376642473, diff(x(t), t) = 54.7872713137340,
y(t) = .908591100591405, diff(y(t), t) = -56.4425175342350]
The velocity is given by the formula
V:=t->sqrt(rhs(p(4)[3])^2+rhs(p(4)[5])^2);
with V(4) equal to 78.66004629802305, agreeing with Table 4.3.5.
===== end