CE 206: Engineering
Computation
p
Sessional
Ordinary Differential Equations
(ODE)
Why study Differential Equations?
Many physical phenomena are best formulated
mathematically in terms of their rate of change.
Motion of a swinging pendulum
Bungee-jumper
Equation
d 2θ g
+ sin θ = 0
2
dt
l
dv
cd 2
=g
g− v
dt
m
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Euler’s Method
The first derivative provides a direct estimate of the
slope at ti:
dy
= f (t i , yi )
dt ti
Euler method uses that
estimate as the increment
ffunction:
φ = f (t i , yi )
yi+1 = yi + f (t i , yi )h
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
MATLAB Code: Euler’s Method
function [t
[t,y]
y] = eulode(tspan,y0,h)
eulode(tspan y0 h)
%
%
%
%
input: dydt = name of the M-file that evaluates the ODE
tspan = [ti, tf], ti and tf limits of independent variable
y0 = initial value of dependent variable
h = step size
output:y = solution vector
if nargin<3,error(‘less than 3 input arguments'),end
ti = tspan(1); tf = tspan(2);
t = (ti:h:tf)'; n = length(t);
% if necessary, add an additional value of t
% so that range goes from t = ti to tf
if t(n)<tf
yi +1 = yi + f ti , yi h
t(n+1) = tf; n = n+1;
end
y = y0*ones(n,1); %preallocate y to improve efficiency
for i = 1:n-1 %implement Euler's method
y(i+1) = y(i) + dydt(t(i),y(i))*(t(i+1)-t(i));
end
(
CE 206: Engg. Computation Sessional
)
Dr. Tanvir Ahmed
Exercise: Euler’s method
Solve the bungee jumper problem using Euler
Euler’ss method (use a
step size of 2 s). Compare the results with the analytical solution
cd 2
dv
=g− v
dt
m
Given g = 9.81
Given,
9 81 m/s2, cd = 0.25
0 25 kg/m,
kg/m m = 68
68.11 kg
60
Analytical solution:
50
v(t ) =
v (m/s)
40
⎛ gcd ⎞
gm
tanh⎜⎜
t ⎟⎟
cd
⎝ m ⎠
30
20
10
0
0
2
4
6
8
10
12
t (sec)
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Heun’s Method
Predictor:
yi0+1 = yi + f (ti , yi )h
f (ti , yi ) + f (ti +1 , yi0+1 )
Corrector: yi +1 = yi +
h
2
Exercise: Modify
y eulode.m to implement
p
Heun’s method.
Solve the same bungee-jumper problem.
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
4th order Runge-Kutta method
1
yi+1 = yi + (k1 + 2k2 + 2k3 + k4 )h
6
where
k1 = f (t i , yi )
⎛
1
1 ⎞
k2 = f ⎜t i + h, yi + k1h ⎟
⎝
2
2 ⎠
⎛
⎞
1
1
k3 = f ⎜t i + h, yi + k2 h ⎟
⎝
⎠
2
2
k4 = f (t i + h,
h yi + k3h)
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
MATLAB Code: Runge-Kutta Method
function [tp,yp] = rk4(tspan,y0,h)
% input:
dydt = name of the M-file that evaluates the ODEs
%
tspan = [ti, tf]; initial and final times
%
y0
0 = i
initial
iti l values
l
of
f d
dependent
d t variables
i bl
%
h = step size
% output: tp = vector of independent variable
%
yp = vector of solution for dependent
p
variables
if nargin<3,error(‘3 input arguments required'), end
ti = tspan(1); tf = tspan(2);
tp = (ti:h:tf)'; n = length(tp);
% if necessary, add an additional value of t
% so that range goes from tp = ti to tf
if tp(n)<tf
t ( )<tf
tp(n+1) = tf; n = n+1;
end
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
MATLAB Code: Runge-Kutta Method
y = y0*ones(n,1); %preallocate y to improve efficiency
for i = 1:n-1 %implement RK method
hh = tp(i+1)-tp(i); tt = tp(i); yy = yp(i);
k1
k2
k3
k4
=
=
=
=
dydt(tt,yy);
dydt(tt + hh/2, yy + 0.5*k1*hh);
dydt(tt + hh/2,
hh/2 yy + 0
0.5*k2*hh);
5*k2*hh);
dydt(tt + hh, yy + k3*hh);
yp(i+1) = yy + hh
hh*(k1
(k1 + 2*k2
2 k2 + 2*k3
2 k3 + k4)/6;
end
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Practice Problem
Solve the following initial value problem over the interval
from t = 0 to 2 where y(0) = 1. Display all your results in
the same graph
dy
= yt 3 − 1.5 y
dt
(a) Analytically
(b) Euler
Euler’ss method with h = 0.5
0 5 and h = 0.25
0 25
(c) Heun’s method with h = 0.5
(d) 4th order RK method with h = 0.5
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
MATLAB’s common ODE solvers
MATLAB’s ode23 function uses second- and third-order
RK functions to solve the ODE and adjust step sizes.
sizes
MATLAB’s ode
ode45 function uses fourth- and fifth-order RK
functions to solve the ODE and adjust step sizes. This is
recommended as the first function to use to solve a
problem.
problem
MATLAB s ode113 function is a multistep solver useful
MATLAB’s
for computationally intensive ODE functions.
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Example of ODE solver: ode45
The functions are generally called in the same way
[t y] = ode45(odefun
[t,
ode45(odefun, tspan,
tspan y0)
y solution array,
y:
y where each column represents
p
one of the
variables
bl and
d each
h row corresponds
d to a time in the
h t vector
odefun: function returning a column vector of the righthand-sides of the ODEs
tspan: time over which to solve the system
If tspan has two entries, the results are reported for those times as
well as several intermediate times based on the steps taken by the
algorithm
If tspan has more than two entries, the results are reported only for
those specific times
y0: vector of initial values
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Solving a system of ODE using ode45
Predator-prey problem
Solve
dy1
= 1.2y
1 2y1 − 0.6y
0 6y1 y2
dt
dy2
= 0.8y
0 8y2 + 0.3y
0 3y1 y2
dt
with y1(0)=2 and y2(0)=1 for 20 seconds
function yp = predprey(t, y)
yp = [1.2*y(1)-0.6*y(1)*y(2);…
-0.8*y(2)+0.3*y(1)*y(2)];
predprey.m
p
p y file
tspan = [0 20];
y0
0 = [2
[2, 1]
1];
[t, y] = ode45(@predprey, tspan, y0);
figure(1); plot(t,y);
figure(2); plot(y(:
plot(y(:,1),y(:,2));
1) y(: 2));
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Example using ode45
Predator-prey problem
Solve
dy1
= 1.2y
1 2y1 − 0.6y
0 6y1 y2
dt
dy2
= 0.8y
0 8y2 + 0.3y
0 3y1 y2
dt
with y1(0)=2 and y2(0)=1 for 20 seconds
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Practice Problem: motion of a pendulum
d 2θ g
+ sin θ = 0
2
dt
l
Analytical solution (assuming θ = sinθ):
⎛ g ⎞
t ⎟⎟
⎝ l ⎠
θ (t ) = θ 0 cos⎜⎜
If l = 2 ft, g = 32.2 ft/s2 and θ0 = π/4, Solve for θ from t = 0 to 1.6 s using
(a) Euler
Euler’ss method with h = 0.05
0 05
(b) Using ode45 function with h = 0.05
(c) Plot your results and compare with the analytical solution
dθ g
+ sin θ = 0
2
dt
l
2
CE 206: Engg. Computation Sessional
dθ
=v
dt
dv
g
= − sin θ
dt
l
@t = 0, θ0 = π/4
@t = 0, v = 0
Dr. Tanvir Ahmed
Boundary value problem: the shooting method
the boundary-value problem is converted into an
equivalent initial-value problem.
⎧ dT
⎪ dx = z
d 2T
+ h′(T∞ − T ) = 0 ⇒ ⎨
2
dz
dx
⎪
= −h′(T∞ − T )
⎩ dx
Generally, the equivalent system will not have
sufficient initial conditions
-A guess is made for any undefined values.
- The guesses are changed until the final solution satisfies all the B.C.
For linear ODEs, only two “shots” are required - the
proper initial condition can be obtained as a linear
i
interpolation
l i off the
h two guesses.
CE 206: Engg. Comp. Sessional
Dr. Tanvir Ahmed
Problem: shooting method
Solve
d 2T
4
4
′
′
+
h
T
−
T
+
σ
T
−
T
=0
(
)
(
)
∞
2
∞
dx
with σσ’=22.7x10
7x10-9 K-3 m-2, L
L=10
10 m,
m hh’=00.05
05 m
m-22, T∞=200
200 K,
K T(0) =
300 K, and T(10) = 400 K.
First - break into two equations:
q
⎧ dT
⎪ dx = z
d 2T
4
4
+ h′(T∞ − T ) + σ ′(T ∞ − T )= 0 ⇒ ⎨
2
d
dx
⎪ dT = −0.05(200 − T ) − 2.7 ×10 −9 (1.6 ×10 9 − T )
⎩ dz
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Problem: shooting method
Code for derivatives:
function dy=dydxn(x,y)
dy=[y(2);…
-0.05*(200-y(1))-2.7e-9*(1.6e9-y(1)^4)];
0 05*(200 (1)) 2 7 9*(1 6 9 (1)^4)]
Code for residual:
function r=res(za)
[x,y]=ode45(@dydxn, [0 10], [300 za]);
r=y(length(x),1)-400;
Code for finding root of residual:
fzero(@res, -50)
Code for solving system:
[x,y]=ode45(@dydxn, [0 10], [300 fzero(@res, -50) ]);
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Practice Problem: shooting method
The b
Th
basic
i diff
differential
i l equation
i off the
h elastic
l i curve ffor a uniformly
if
l lloaded
d d
beam is given as
d 2 y wLx wx 2
EI 2 =
−
dx
2
2
If E = 30000 ksi, I = 800 in4, w = 1 kip/ft,
p/ L = 10 ft, solve for the deflection
of the beam using the shooting method and compare the numerical
results with the analytical solution
L 3 wx 4 wL
L3 x
wLx
y=
−
−
12 EI 24 EI 24 EI
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed
Numerical solution: finite difference
finite differences are substituted for the derivatives
in the original equation
d 2T
+ h′(T∞ − T ) = 0
2
dx
d 2T Ti−1 − 2Ti + Ti+1
=
2
dx
Δx 2
Ti−1 − 2Ti + Ti+1
+ h′(T∞ − Ti ) = 0
2
Δx
−Ti−1 + (2 + h′Δx 2 )Ti − Ti+1 = h′Δx 2T∞
CE 206: Engg. Comp. Sessional
Dr. Tanvir Ahmed
Numerical solution: finite difference
The linear
Th
li
diff
differential
ti l equation
ti iis ttransformed
f
d iinto
t a
set of simultaneous algebraic equations.
Since T0 and Tn are known (Drichlet boundary
conditions), they will be on the right-hand-side of the
linear algebra system
⎡2 + h′Δx 2
⎢ −1
⎢
⎢
⎢⎣
−11
2 + h′Δx 2
O
CE 206: Engg. Comp. Sessional
−1
O
−1
⎤⎧ T ⎫ ⎧h′Δx 2T∞ + T0 ⎫
⎥⎪ T1 ⎪ ⎪⎪ h′Δx 2T ⎪⎪
∞
⎥⎨ 2 ⎬ = ⎨
⎬
O ⎥⎪ M ⎪ ⎪
M
⎪
2 + h′Δx 2 ⎥⎦⎩Tn−1 ⎭ ⎪⎩h′Δx 2T∞ + Tn ⎪⎭
Dr. Tanvir Ahmed
Numerical solution: finite difference
If there
h
iis a d
derivative
i i b
boundary
d
condition
di i
(Neumann B.C.), the centered difference equation is
solved at the point and rewriting the system
equation accordingly.
For a Neumann condition at T0 point:
⎛ dT ⎞
dT
T1 − T−1
=
⇒ T−1 = T1 − 2Δx⎜
⎟
2Δx
dx 0
dx
⎝
0⎠
−T−1 + (2 + h′Δx 2 )T0 − T1 = h′Δx 2T∞
⎡
⎛ dT ⎞⎤
2
2
−⎢T1 − 2Δx⎜
⎟⎥ + (2 + h′Δx )T0 − T1 = h′Δx T∞
⎝ dx 0 ⎠⎦
⎣
⎛ dT ⎞
2
2
(2 + h′Δx )T0 − 2T1 = h′Δx T∞ − 2Δx⎜⎝ dx ⎟⎠
0
CE 206: Engg. Comp. Sessional
Dr. Tanvir Ahmed
Practice Problem: Finite difference
Solve the nondimensionalized ODE using finite difference
methods that describe the temperature distribution in a
circular rod with internal heat source S
d 2T 1 dT
+
+S =0
2
dr
r dr
Over the range 0 ≤ r ≤ 1, with boundary conditions
T (r = 1) = 1
dT
dr
=0
r =0
For S = 1, 10 and 20 K/m2 plot the temperature versus radius
CE 206: Engg. Computation Sessional
Dr. Tanvir Ahmed

Download Report