Transcript Day 23 slides

Scientific Computing
Systems Of
Ordinary Differential Equations
One-D System
• The ODE +IV
dx
 f t, x 
dt
x(a)= x0 = c
has one independent variable – t – and one
dependent
 variable – x. This is a onedimensional ODE system.
Higher Dimensional Systems
•
dx
 f t , x(t ), y (t ) 
dt
dy
 g t , x(t ), y (t ) 
dt
x(a)  c
y (a)  d
is a two-dimensional ODE system. A solution
is a vector function (x(t), y(t)) that satisfies
this system.
Higher Dimensional Systems
• Example:
dx
 t2  x
dt
dy
 y2  y  t
dt
x(0)  1
y (0)  0
is actually solvable as two separate 1-d
systems. A system that can be separated like
this is called an uncoupled system.
Higher Dimensional Systems
• Example:
dx
 t2  x  y
dt
dy
 y2  x  t
dt
x(0)  1
y (0)  0
is not solvable as two separate 1-d systems. A
system that cannot be separated is called a
coupled system.
Higher Dimensional Systems
• In General we have n-dimensional system:
dx1
 f1 (t , x1 (t ), x2 (t ), , xn (t ))
dt
dx2
 f 2 (t , x1 (t ), x2 (t ), , xn (t ))
dt

dxn
 f n (t , x1 (t ), x2 (t ), , xn (t ))
dt
x1 ( a )  c1 ,  , xn (a )  cn
Higher Dimensional Systems
• In vector form, if:
then,
 x1 
 f1 
 c1 
 
 
 
x
f
c2 
2
2



X
, F
, C
 
  

 
 
 
 xn 
 f n 
cn 
d
X (t )  F(t , X (t )),
dt
X ( a)  C
Note the similar form to the 1-d system case.
Lotka VolterraSystems
• The Lotka-Volterra system models the
interactions between two species in an
ecosystem, a predator and a prey (foxes and
rabbits). Let r(t) and f(t) represent the number
of rabbits and foxes that are alive at time t.
The Lotka-Volterra model is:
dr
 ar  b rf
dt
df
 eb rf  cf
dt
Lotka VolterraSystems
• The Lotka-Volterra model is:
where
•
•
•
•
dr
 ar  b rf
dt
df
 eb rf  cf
dt
a = natural growth rate of rabbits (absent predators)
c = natural death rate of foxes in the absence of food (rabbits)
b = death rate per encounter of rabbits with foxes
e = efficiency of turning gobbled rabbits into foxes.
Solving Systems
• To solve
d
X (t )  F(t , X (t )),
dt
X ( a)  C
we can use the various algorithms already
discussed, only casting them in vector form.
Vector Taylor’s Series:
h2
hk (k )
X (t0  h) X (t0 )   h X ' (t0 ) 
X ' ' (t 0 )   
X (t 0 )  
2
k!
Euler Method for Systems
X ti 1   X ti   hFti , X (ti )
 X ti   hFti , x1 ti , x2 (ti ),, xn (ti )
Heun’s Method for Systems
1
1

X i 1  X i   K1  K 2 h
2 
2
K1  Fti , X i 
K 2  Fti  h, Xi  K1h
Runge-Kutta Order 4 Method
Xi 1  Xi  φi h
φi  (K1  2K 2  2K3  K 4 ) / 6
K1  F(ti , X i )
K 2  F(ti  h / 2, Xi  K1h / 2)
K3  F(ti  h / 2, Xi  K 2h / 2)
K 4  F(ti  h, Xi  K 3h)
Converting Higher Order ODES
to Systems
• A general, nth-order ODE
( n)
( n 1)

x

f
(
t
,
x
,
x
'
,
x
'
'
,

,
x
)


( n 1)
x
(
t
)


,
x
'
(
t
)


,

,
x
(t0 )   n 1

0
0
1
 0
can be converted to a system of ODE’s
 x1  x
x1 ( t 0 )  α 0
 x1 '  x2 ,


x2 ( t 0 )  α1
x

x
'
 2
 x2 '  x3 ,


let  x3  x' '
  x3 '  x4 ,
x3 ( t 0 )  α 2
 
 


 x  x ( n 1) 
 xn '  f ( t, x1, x 2 ,, x n ), xn ( t 0 )  α n 1
 n
Example:Bungee Jumper
Coupled ODEs
Bungee Jumper – Coupled ODEs
• Forces:
– mg (gravity, g = acceleration due to gravity)
– cd v2 (drag force, cd = drag coefficient, v = velocity)
(need to always retard v, so if falling (v>0) need force neg,
if rising (v<0) need force pos to reduce dv/dt)
(use sign(v) for drag force)
– k (x-L) (spring force, x = distance measured down from
platform, L = rest length of cord)
– γv
(damping force, γ is damping coefficient of cord)
Bungee Jumper – Coupled ODEs
dv
k =  = 0 if
2
m
 mg  sign(v)c d v  k(x  L)  v x  L
dt
• Dividing by m:
cd 2 k
dv

 g  sign( v ) v  ( x  L)  v
dt
m
m
m
k =  = 0 if
xL
Bungee Jumper – Coupled ODEs
• Convert to a system of ODEs for x and v
 dx
 dt  v

 dv  g  sign( v ) c d v 2  k ( x  L)   v
 dt
m
m
m
k =  = 0 if
xL
Matlab Vector Runge4
Change code how??
function [ t x ] = runge(f, a, c, h, n )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
x = zeros(n+1,1);
x(1) = c; t = [a:h:n*h + a];
for i = 1:n
k1 = f(t(i),x(i));
k2 = f(t(i) + h/2, x(i) + k1*h/2);
k3 = f(t(i) + h/2, x(i) + k2*h/2);
k4 = f(t(i) + h, x(i) + k3*h);
phi = (k1 + 2*k2 + 2*k3 + k4)/6;
x(i+1) = x(i) + phi*h;
end
t=t';