E7_L27_ODE2_Euler_mod_Euler_F08

E7_L27_ODE2_Euler_mod_Euler_F08 - 1 E7: INTRODUCTION TO...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

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 (predictor-corrector) 3. Examples Copyright 2007, Horowitz, Packard. This work is licensed under the Creative Commons Attribution-Share 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/by-sa/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 first-order 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 Runge-Kutta (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 right-hand-side 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 (predictor-corrector) – 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 pseudo-code 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 row-wise (continues) E7 L27 Euler pseudo-code 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:M-1 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 second-order 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 second-order 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 pseudo-code 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 pseudo-code …. 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:M-1 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 second-order approximation update % store results row-wise (continues) E7 L27 SetSet-up 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 first-order 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 second-order term of the Taylor series expansion • Runge-Kutta methods (high accuracy) – utilize up to the fourth-order 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.

Ask a homework question - tutors are online