Scribd
Upload
392 views23 pages

Ce206 Ode

This document discusses numerical methods for solving ordinary differential equations (ODEs) including Euler's method, Heun's method, the fourth-order Runge-Kutta method, and MATLAB implementations. It also covers solving boundary value problems using the shooting method and finite difference methods. Some key topics covered are: 1. Euler's method, Heun's method, and the 4th order Runge-Kutta method for solving initial value problems. 2. MATLAB code examples for implementing these methods. 3. Using MATLAB's ode45 and ode113 solvers to solve problems. 4. The shooting method for converting boundary value problems to initial value problems. 5. Finite difference methods

Uploaded by

Borhan Shaikat
Copyright
© Attribution Non-Commercial (BY-NC)
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
392 views23 pages

Ce206 Ode

This document discusses numerical methods for solving ordinary differential equations (ODEs) including Euler's method, Heun's method, the fourth-order Runge-Kutta method, and MATLAB implementations. It also covers solving boundary value problems using the shooting method and finite difference methods. Some key topics covered are: 1. Euler's method, Heun's method, and the 4th order Runge-Kutta method for solving initial value problems. 2. MATLAB code examples for implementing these methods. 3. Using MATLAB's ode45 and ode113 solvers to solve problems. 4. The shooting method for converting boundary value problems to initial value problems. 5. Finite difference methods

Uploaded by

Borhan Shaikat
Copyright
© Attribution Non-Commercial (BY-NC)
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as PDF, TXT or read online on Scribd
You are on page 1/ 23

CE 206: Engineering Computation Sessional p 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 v g dt m

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

Eulers 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 function: f

= f (t i , yi )
yi+1 = yi + f (t i , yi )h

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

MATLAB Code: Eulers Method


function [t y] = eulode(tspan y0 h) [t,y] 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: Eulers method


Solve the bungee jumper problem using Euler s method (use a Eulers step size of 2 s). Compare the results with the analytical solution

cd 2 dv =g v dt m
60

Given, Given g = 9.81 m/s2, cd = 0.25 kg/m, m = 68 1 kg 9 81 0 25 kg/m 68.1

50

Analytical solution:
v(t ) = gcd gm tanh m t cd

40

v (m/s)

30

20

10

10

12

t (sec)
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Heuns 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 eulode.m to implement Heuns method. y p 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 yi + k3h) h,

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 = i iti l values of d 0 initial l f dependent variables d t i bl % h = step size % output: tp = vector of independent variable % yp = vector of solution for dependent variables p

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 t ( )<tf tp(n)<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 yy + 0 5*k2*hh); hh/2, 0.5*k2*hh); dydt(tt + hh, yy + k3*hh);

yp(i+1) = yy + hh*(k1 + 2*k2 + 2*k3 + k4)/6; hh (k1 2 k2 2 k3 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) Eulers method with h = 0.5 and h = 0.25 Euler s 05 0 25 (c) Heuns method with h = 0.5 (d) 4th order RK method with h = 0.5

CE 206: Engg. Computation Sessional

Dr. Tanvir Ahmed

MATLABs common ODE solvers


MATLABs ode23 function uses second- and third-order RK functions to solve the ODE and adjust step sizes sizes. MATLABs ode45 function uses fourth- and fifth-order RK ode 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 MATLABs ode113 function is a multistep solver useful 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, [t y] = ode45(odefun tspan y0) ode45(odefun, tspan,
y y: solution array, where each column represents one of the y p variables and each row corresponds to a time in the t vector bl d h d h 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 2y1 0 6y1 y2 1.2y 0.6y dt dy2 = 0 8y2 + 0 3y1 y2 0.8y 0.3y 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)];

p p y predprey.m file

tspan = [0 20]; y0 = [2 1] 0 [2, 1]; [t, y] = ode45(@predprey, tspan, y0); figure(1); plot(t,y); figure(2); plot(y(: 1) y(: 2)); plot(y(:,1),y(:,2));
CE 206: Engg. Computation Sessional Dr. Tanvir Ahmed

Example using ode45


Predator-prey problem Solve
dy1 = 1 2y1 0 6y1 y2 1.2y 0.6y dt dy2 = 0 8y2 + 0 3y1 y2 0.8y 0.3y 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):

(t ) = 0 cos

g t l

If l = 2 ft, g = 32.2 ft/s2 and 0 = /4, Solve for from t = 0 to 1.6 s using (a) Eulers method with h = 0 05 Euler s 0.05 (b) Using ode45 function with h = 0.05 (c) Plot your results and compare with the analytical solution
d =v dt dv g = sin dt l

d g + sin = 0 2 dt l
2

@t = 0, 0 = /4 @t = 0, v = 0
Dr. Tanvir Ahmed

CE 206: Engg. Computation Sessional

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 interpolation of the two guesses. i l i f h
CE206:Engg.Comp.Sessional Dr.Tanvir Ahmed

Problem: shooting method


Solve
d 2T + h(T T ) + ( 4 T 4 )= 0 T 2 dx

with =2 7x10-9 K-3 m-2, L=10 m h=0 05 m-2 T=200 K T(0) = 2.7x10 L 10 m, h 0.05 m 2, 200 K, 300 K, and T(10) = 400 K.

First - break into two equations: q


dT dx = z d 2T 4 4 + h(T T ) + ( T )= 0 T 2 d dx dT = 0.05(200 T ) 2.7 10 9 ( 10 9 T ) 1.6 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 basic differential equation of the elastic curve f a uniformly l d d Th b i diff i l i f h l i for if l loaded 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, L = 10 ft, solve for the deflection p/ of the beam using the shooting method and compare the numerical results with the analytical solution

L L wLx 3 wx 4 wL3 x 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 Ti1 2Ti + Ti+1 = 2 dx x 2 Ti1 2Ti + Ti+1 + h(T Ti ) = 0 2 x Ti1 + (2 + hx 2 ) i Ti+1 = hx 2T T


CE 206: Engg. Comp. Sessional Dr. Tanvir Ahmed

Numerical solution: finite difference


The li Th linear diff differential equation i t ti l ti is transformed i t a f d into 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 + hx 2 1 1 1 2 + hx 2 O 1 O 1 T hx 2T + T0 T1 hx 2T 2 = O M M 2 + hx 2 Tn1 hx 2T + Tn

CE 206: Engg. Comp. Sessional

Dr. Tanvir Ahmed

Numerical solution: finite difference


If there i a d i i b h is derivative boundary condition d 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 T1 = T1 = T1 2x 2x dx 0 dx 0 T1 + (2 + hx 2 ) 0 T1 = hx 2T T dT 2 2 T1 2x T + (2 + hx ) 0 T1 = hx T dx 0 dT 2 2 (2 + hx )T0 2T1 = hx T 2x 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

You might also like