This preview shows page 1. Sign up to view the full content.
Unformatted text preview: 1 E7: INTRODUCTION TO COMPUTER PROGRAMMING FOR SCIENTISTS AND PROGRAMMING FOR SCIENTISTS AND ENGINEERS ENGINEERS
Lecture Outline Review of Ordinary Differential Equations (ODEs) 2 Many problems in science and engineering lead to problems in science and engineering lead to Ordinary Differential Equations (ODEs) of the form 1. 2. Review of ordinary differential equations Numerical integration algorithms: – Euler – Modified Euler (predictorcorrector) 3. Examples
Copyright 2007, Horowitz, Packard. This work is licensed under the Creative Commons AttributionShare 2007 Horowitz Packard This work is licensed under the Creative Commons Attribution Alike License. To view a copy of this license, visit http://creativecommons.org/licenses/bysa/2.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA. E7 L27 dy (t ) = f (t , y ) dt
where: • t is the independent scalar variable (often time) • y is the dependent variable, which can be a vector • f(t,y) is a known function of its arguments
E7 L27 Solution of ODEs (Initial Value Problem) • Given the ODE the ODE 3 What is the difference between ode45 and quad? Quad is used to solve integrals. For example: is used to solve integrals For example: 4 • and an initial condition • find the function dy (t ) = f (t , y ) dt y (t0 ) = y0
(slight abuse of notation) dy (t ) = f (t ) dt
Only a function of time t This can be solved by quadratures: y (t ) y (t0 ) = y0
• which satisfies: satisfies: dy (t ) = f (t , y (t )) dt for t ≥ t0
E7 L27 y (t ) = ∫ f (τ )dτ + y (to )
to
E7 L27 t What is the difference between ode45 and quad? Quad is used to obtain definite integrals numerically is used to obtain definite integrals numerically (e.g. quadratures): 5 What is the difference between ode45 and quad? ode45 is used to solve ODEs: is used to solve ODEs: 6 q = ∫ f (t )dt
a
Quad syntax: b dy (t ) = f (t , y ) dt
a function of time t AND y This cannot solved by quad, unless we can use cannot separation of variables, e.g. q = quad(@f,a,b) quad(@
where f is only a function of time.
E7 L27 f ( t ) = g ( t ) h( y )
E7 L27 Higher order ODEs Remember the basic ODE: the basic ODE: 7 Higher order ODEs System of n firstorder coupled ODEs : ODE 8 y = f (t , y )
y
where: and y = f (t , y )
y can be vectors!
then f will also be a vector • t is the independent scalar variable (often time) • y is the dependent variable • f(t,y) is a known function of its arguments
E7 L24 y2 = f 2 ( t , y1 , y2 ,..., yn ) yn = f n ( t , y1 , y2 ,..., yn )
E7 L24 y1 = f1 ( t , y1 , y2 ,..., yn ) Euler’s Euler’s disk (ME170) – courtesy Prof. O’Reilly
Define: 9 Second order pendulum example (Palm 8.6–1) A physics 7A and ME104 classical problem physics 7A and ME104 classical problem 10 y1 = α ,
⎡ ⎢ ⎢ ⎣ y2 = ω1 ,
⎡ ⎢ ⎢ ⎣ h y3 = ω2 ,
y2 y4 = ω3 ,
i⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ y1 ⎢ ⎥ ⎢ 1 1 y (2y − y tan(y )) − 2 y 2 sin(2y ) + sin(y ) ⎢y ⎥ ⎢ 3 4 1 1 1 4 2 ⎢ 54 52 d⎢ 2⎥ ⎢ ⎥ = ⎢ (1− 5 cos (y1)) ⎢y ⎥ ⎢ 0 dt ⎢ 3 ⎥ ⎢ ⎤ ⎥ ⎥ ⎦ θ+ g sin(θ ) = 0 L
second order nonlinear ODE ODE L θ y4 y2(y4 tan(y1) − 2y3) g http://www.youtube.com/watch?v=9SUA7C0cBd8 http://www.fpp.edu/~milanb/euler/ rotation θ angular rotation θ angular acceleration
E7 L24 E7 L24 Second order pendulum example (Palm 8.6–1) 11 Numerical solution of ODEs using Matlab Matlab has a variety of ODEs numerical integrators. has variety of ODEs numerical integrators The two that we will be using here are: 12 d y1 = y2 dt d g y2 = − sin( y1 ) dt L
equivalent descriptions! L θ • ode45: an explicit RungeKutta (RK) algorithm.
In general, the best function to apply as a "first try" for most problems. g • ode23: also an explicit RK algorithm.
More efficient at crude tolerances than ode45 in the efficient crude tolerances than ode45 in the presence of moderate stiffness. Both functions have the same syntax. g θ + sin(θ ) = 0 L y1 = θ angular rotation y2 = θ angular velocity
E7 L24 E7 L27 ode45 basic syntax (for more options see help) 13 ode45 basic syntax (for more options see help) 14 [T,Y] = ode45(@odefun,tspan,y0); ode45(@odefun
input arguments: [T,Y] = ode45(@odefun,tspan,y0); ode45(@odefun
output arguments: • odefun: function that evaluates righthandside
of ODE • T: column vector of time points ti • Y : solution array: Y(n) = y(T(n))’ ? y = f (t , y ) • tspan: normally, tspan = [t0, tfinal] • y0: initial condition y0 = y(t0)
E7 L24 y = f (t , y )
and y can be vectors! (more later)
E7 L24 y Odefun – ydot_pend1 for pendulum example 15 Using ode45 – use anonymous function mgl = 9.81; ydot_pend = @(t,y) ydot_pend1(y,mgl);
handle to anonymous function with input arguments t and y 16 y1 = y2
g y2 = − sin( y1 ) L
⎡y ⎤ y = ⎢ 1⎥ ⎣ y2 ⎦ y1 = θ
y2 = θ θ L g
2 ×1 arrays (column vectors) ydot_pend1.m ⎡y ⎤ y = ⎢ 1⎥ ⎣ y2 ⎦ [T,Y]=ode45(ydot_pend,tspan,y0);
ydot_pend1.m function ydot = ydot_pend1(y,mgl) % mgl = g/L ydot = [ y(2); mgl * sin(y(1))];
E7 L24 function ydot = ydot_pend1(y,mgl) % mgl = g/L ydot = [ y(2); mgl * sin(y(1))];
E7 L24 Using ode45 – setting y0 and tspan y0 tspan mgl = 9.81; ydot_pend = @(t,y) ydot_pend1(y,mgl); y0=[pi/4;0]; tspan=[0;5];
⎡π / 4 ⎤ y (0) = ⎢ ⎥ ⎣0⎦ 17 Using ode45 on pendulum example
[ T , Y ] = ode45(ydot pend, tspan, y0); ode45(ydot_pend tspan y0); subplot(2,1,1) plot(T,Y(:,1)) xlabel('time  seconds') ylabel('\theta  rad') 18 ⎡ y1 (0)⎤ ⎡θ (0)⎤ ⎢ ⎥⎢ ⎥ Y (:,1) = ⎢ ⎥=⎢ ⎥ ⎢ y1 (5)⎥ ⎢θ (5) ⎥ ⎣ ⎦⎣ ⎦ [T,Y]=ode45(ydot_pend,tspan,y0);
ydot_pend1.m function ydot = ydot_pend1(y,mgl) % mgl = g/L ydot = [ y(2); mgl * sin(y(1))];
E7 L24 ⎡ y2 (0)⎤ ⎡θ (0)⎤ subplot(2,1,2) ⎥ ⎥=⎢ Y (:,2) = ⎢ plot(T,Y(:,2)) ⎢ ⎥ ⎢ ⎥ ⎢ y2 (5) ⎥ ⎢θ (5) ⎥ xlabel xlabel('time  seconds') ⎣ ⎦⎣ ⎦ ylabel('d/dt \theta  rad/sec') E7 L24 Using ode45 on pendulum example [ T , Y ] = ode45(ydot_pend, tspan, y0); ode45( tspan y0);
subplot(2,1,1) plot(T,Y(:,1)) xlabel('time  seconds') ylabel('\theta  rad')
θ  rad 19 Numerical integration of ODEs Today we will explain we will explain • Euler’s method
– Approximates a Taylor expansion up to its first order Taylor expansion up to its first order 20 1 0.5 0 0.5 1 0 y1 = θ
1 2 3 time  seconds 4 5 subplot(2,1,2) plot(T,Y(:,2)) xlabel('time  seconds') ylabel('d/dt \theta  rad/sec') • Modified Euler’s method (predictorcorrector)
– Approximates a Taylor expansion up to its second order Taylor expansion up to its second order d/dt θ  rad/sec y1 = θ
y2 = θ 4 θ L g 2 0 2 4 0 y2 = θ
1 2 3 time  seconds 4 5 E7 L24 E7 L27 Euler’s method Consider the ODE the ODE 21 Taylor’s series expansion 22 y (t ) = dy (t ) = f (t , y ) dt y (t1 ) = y (t0 + h)
can be obtained from the Taylor series expansion: be obtained from the Taylor series expansion: and the initial condition y (t0 ) = y0 y (t0 + h ) = y (t0 ) + y (t0 ) ⋅ h +
where where 1 y (t0 ) ⋅ h 2 + … 2 We want to obtain an estimate of want to obtain an estimate of y (t1 ) t1 = t0 + h
and h is the integration step size. integration
E7 L27 E7 L27 Taylor’s series expansion 23 Euler’s method 24 y (t1 ) = y (t0 + h)
can be obtained from the Taylor series expansion: be obtained from the Taylor series expansion: y (t1 ) = y (t0 + h)
can be approximated using only the first order term be app te of the Taylor series expansion: y (t0 + h ) = y (t0 ) + y (t0 ) ⋅ h + ... + 1 [n] y (t0 ) ⋅ h n + … n! y (t0 + h) ≈ y (t0 ) + y (t0 ) ⋅ h
if h is sufficiently small. h
E7 L27 E7 L27 Euler’s method 25 Euler’s method The first step of Euler’s method uses the approximation: 26 y (t0 + h) ≈ y (t0 ) + y (t0 ) ⋅ h
we now substitute the ODE y (t0 ) = f (t0 , y (t0 )) y (t0 + h) ≈ y (t0 ) + f (t0 , y (t0 )) ⋅ h
Thus: in the approximation to obtain: y (t0 + h) ≈ y (t0 ) + f (t0 , y (t0 )) ⋅ h
E7 L27 t1 = t0 + h
y E (t1 ) = y (t0 ) + f (t0 , y (t0 )) ⋅ h
E7 L27 Euler’s method Notice that we are using the superscript E to denote th th the Euler approximate solution of the ODE 27 Euler’s method Notice that we are using the superscript E to denote th th the Euler approximate solution of the ODE and so on…. 28 t1 = t0 + h
y E (t1 ) = y (t0 ) + f (t0 , y (t0 )) ⋅ h
The formula is repeated in the next step: tk +1 = tk + h
y E (tk +1 ) = y E (tk ) + f (tk , y E (tk )) ⋅ h
Initial condition:
E7 L27 t2 = t1 + h
y E (t2 ) = y E (t1 ) + f (t1 , y E (t1 )) ⋅ h y E (t0 ) = y (t0 )
E7 L27 Graphical Graphical explanation of Euler’s method 29 Graphical explanation of Euler’s method 30 y (t ) y (t ) = f (t , y ) y (t ) y E (t2 ) y (t1 ) y E (t1 ) y (t0 ) h y (t ) = f (t , y ) y (t1 ) y E (t1 ) f (t1 , y E (t1 )) ⋅ h y (t0 )
h
t0 t1 f (t0 , y (t0 )) ⋅ h h
t2 t
E7 L26 h
t1 t2 t
E7 L27 t0 Euler’s method: • Start with: with: 31 Euler pseudocode pseudofunction [ T , Y ] = odeeuler( fh, tspan, y0, h) odeeuler( fh tspan y0
% uses same first 3 input arguments as ode45 32 y E (t0 ) = y (t0 ) • Select a sufficiently small integration step size sufficiently small integration step size • Iterate for h % plus integration step size h % fh is function handle to f(t , y) (ydot = f(t,y))
n = length(y0);
% assume yo is a column vector and determine % the order n (dimension of y0) k ≥0 tk +1 = tk + h
y E (tk +1 ) = y E (tk ) + f (tk , y E (tk )) ⋅ h
• Until T = (tspan(1):h:tspan(end))’; % define T (time) column vector M = length(T); % determine number of iterations
% use same convention as ode45 to store results for Y, tk ≥ t final Y = zeros(M,n); zeros(M Y(1,:) = y0’;
E7 L27 % store results rowwise (continues)
E7 L27 Euler pseudocode pseudofunction [ T , Y ] = odeeuler( fh, tspan, y0, h) odeeuler( fh tspan y0
% uses same first 3 input arguments as ode45 % plus integration step size h 33 Euler integration of pendulum example 34 y1 = y2 g y2 = − sin( y1 ) L y1 = θ
y2 = θ θ L % fh is function handle to f(t , y) (ydot = f(t,y)) (from previous) Anonymous function handle ydot_pend(t,y): mgl = 9.81; ydot_pend = @(t,y) ydot_pend1(y,mgl);
ydot_pend1.m g % perform Euler’s numerical integration algorithm for k=1:M1 Y(k+1,:) = Y(k,:) + h* fh(T(k), Y(k,:)’)’; Y(k,:) fh(T(k), Y(k,:) end
E7 L27 function ydot = ydot_pend1(y,mgl) % mgl = g/L ydot = [ y(2); mgl * sin(y(1))];
E7 E7 L27 Using ode45 on pendulum example
[ T , Y ] = ode45(ydot pend, [0,5], [pi/4,0]); t_ [0 35 Euler integration of pendulum example h =.05;
[ T , Y ] = odeeuler(ydot_pend, [0,5], [pi/4;0],h); 36 1 0.5
θ  rad 0 0.5 1 0 θ  rad y1 = θ
1 ode45 2 3 time  seconds 4 5 serious deterioration with increased integration set size h L g 5 0 y1 = θ h=0.05 h=0.1 ode45 5 10 0 10 d/dt θ  rad/sec 5 0 5 10 0 1 2 3 time  seconds 4 5 4 d/dt θ  rad/sec 2 0 2 4 0 y2 = θ
y1 = θ y2 = θ
1 ode45 ode45 h=0.1
2 3 time  seconds h=0.05 y1 = θ
y2 = θ θ L g ode45
1 2 3 time  seconds 4 5 y2 = θ
E7 L27 θ 4 5 E7 L27 Improving on Euler’s algorithm • Euler’s method is based on using the first term of method is based on using the first term of the Taylor series expansion: 37 Modified Euler’s method (second order): 1) Do a function call and an Euler prediction 38 y (t0 + h) ≈ y (t0 ) + y (t0 ) ⋅ h
• Why not use the first and second terms of the Taylor series expansion? f1 = f (tk , y E (tk ))
E yP (tk +1 ) = y E (tk ) + f1 ⋅ h 2) Do another function call and a secondorder update y (t0 + h) ≈ y (t0 ) + y (t0 ) ⋅ h + 1 y (t0 ) ⋅ h 2 2
E7 L27 y E (tk +1 ) = y E (tk ) + h E ⎡ f (tk +1 , yP (tk +1 )) + f1 ⎤ ⎦ 2⎣
E7 L27 Modified Euler’s method (second order): 39 Modified Euler’s method (second order): 40 y E (tk +1 ) = y E (tk ) +
E f (tk +1 , yP (tk +1 )) ⋅ h h E ⎡ f (tk +1 , yP (tk +1 )) + f1 ⎦ ⎤ 2⎣ E f (tk +1 , yP (tk +1 )) ⋅ h E yP (tk +1 ) = y E (tk ) + f1 ⋅ h y E ( tk ) f1 ⋅ h = f (tk , y E (tk )) ⋅ h
E7 L27 h E ⎡ f (tk +1 , yP (tk +1 )) + f1 ⎤ ⎣ ⎦ 2 y E ( tk ) f1 ⋅ h = f (tk , y E (tk )) ⋅ h
E7 L27 Modified Euler’s method (second order): 1) Do a function call and an Euler prediction Do function call and an Euler prediction 41 Modified Euler’s method (second order): • Start with: with: 42 y E (t0 ) = y (t0 ) f1 = f (tk , y E (tk ))
E yP (tk +1 ) = y E (tk ) + f1 ⋅ h • Select the integration step size • Iterate for h k ≥0 tk +1 = tk + h f1 = f (tk , y E (tk ))
E yP (tk +1 ) = y E (tk ) + f1 ⋅ h 2) Do another function call and a secondorder update y E (tk +1 ) = y E (tk ) + h E ⎡ f (tk +1 , yP (tk +1 )) + f1 ⎤ ⎣ ⎦ 2
E7 L27 y E (tk +1 ) = y E (tk ) +
• Until h E ⎡ f (tk +1 , yP (tk +1 )) + f1 ⎤ ⎣ ⎦ 2
E7 L27 tn ≥ t final Modified Euler pseudocode pseudofunction [ T , Y ] = odemodeuler( fh, tspan, y0, h) odemodeuler( fh tspan y0
% uses same first 3 input arguments as ode45 % plus integration step size h 43 Modified Euler pseudocode …. continuation pseudoh2 = h/2 h/2 % define h/2 define h/2 44 % fh is function handle to f(t , y) (ydot = f(t,y))
n = length(y0); length(y0);
% assume yo is a column vector and determine assume yo is column vector and determine % the order n (dimension of y0) % perform Modified Euler’s numerical integration algorithm perform Modified Euler numerical integration algorithm for k=1:M1 f1= fh(T(k), Y(k,:)); yp= Y(k,:)’ + h*f1;
% evaluate first function call % do an Euler prediction T = (tspan(1):h:tspan(end))’; % define T (time) column vector M = length(T); % determine number of iterations
% use same convention as ode45 to store results for Y, Y = zeros(M,n); Y(1,:) = y0’; % Do a secondorder approximation update
% store results rowwise (continues)
E7 L27 SetSetup part is exactly like Euler’s algorithm Y(k+1,:) = Y(k,:) + h2 * (fh(T(k+1), yp) + f1 )'; Y(k h2 (fh(T(k+1) yp) f1 end
E7 L27 Comparison using the pendulum example Numerical solution obtained with ode45: 45 Comparison using the pendulum example Numerical solution obtained with ode45, Euler, and modified Euler. 46 1 0.5
θ  rad 0 0.5 1 0 θ  rad y1 = θ
1 h = 0.05 ode45 modified Euler Euler almost as good as ode45
4 5 3 2 1 0 1 2 0 10 d/dt θ  rad/sec 5 y1 = θ
modified Euler
1 Euler 2 3 time  seconds 2 3 time  seconds 4 5 4 d/dt θ  rad/sec 2 0 2 4 0 y2 = θ
y1 = θ y2 = θ
modified Euler
1 Euler y1 = θ y2 = θ θ L g ode45
1 2 3 time  seconds 4 5 y2 = θ
E7 L27 θ L g 0 5 0 2 3 time  seconds 4 5 E7 L27 Comparison using the pendulum example Numerical solution obtained with ode45, Euler, and modified Euler. 47 Summary What did we learn today? did we learn today? Numerical solutions of ODEs can be obtained using a variety of methods: • Euler's method (low accuracy) – utilizes the firstorder term of the Taylor series expansion 48 h = 0.1 modified Euler Euler is significantly more accurate than Euler L g
θ  rad 5 0 y1 = θ
modified Euler
1 5 Euler
4 5 10 0 5 d/dt θ  rad/sec 0 2 3 time  seconds y2 = θ
modified Euler
1 • Modified Euler's method (better accuracy) – utilizes up to the secondorder term of the Taylor series expansion • RungeKutta methods (high accuracy) – utilize up to the fourthorder term of the Taylor series expansion y1 = θ y2 = θ θ 5 Euler
4 5 10 0 2 3 time  seconds E7 L27 E7 L27 ...
View
Full
Document
This note was uploaded on 11/01/2009 for the course ENGLISH 7 taught by Professor Sengupta during the Spring '09 term at University of California, Berkeley.
 Spring '09
 Sengupta
 Ode

Click to edit the document details