1y ago

37 Views

4 Downloads

4.71 MB

60 Pages

Transcription

SECTION 9: ORDINARYDIFFERENTIAL EQUATIONSMAE 4020/5020 – Numerical Methods with MATLAB

2K. WebbIntroductionMAE 4020/5020

Ordinary Differential Equations3 Differential equations can be categorized as eitherordinary or partial differential equations Ordinary differential equations (ODE’s) – functions of asingle independent variable Partial differential equations (PDE’s) – functions of two ormore independent variablesWe’ll focus on ordinary differential equations onlyNote that we are not making any assumption oflinearity here K. WebbAll techniques we’ll look at apply equally to linear ornonlinear ODE’sMAE 4020/5020

Differential Equation Order4 The order of a differential equation is the highestderivative it containsFirst‐order ODE’s contain only first derivatives Second‐order ODE’s include second derivatives (possiblyfirst, as well), and so on Any‐ order ODE can be reduced to a system offirst‐order ODE’s Solution requires knowledge ofconditionsinitial or boundaryWe’ll focus on techniques to solve first‐order ODE’s K. WebbCan be applied to systems of first‐order ODE’s representinghigher‐order ODE’sMAE 4020/5020

Initial‐Value vs. Boundary‐Value Problems5 To solve an‐order ODE (or a system of first‐order ODE’s), known conditions are required Initial‐valueproblems – all conditions are specified atthe same value of the independent variable (typically,ator) Boundary‐valueproblems – conditions specified atdifferent values of the independent variable In this course, we’ll focus exclusively on initial‐valueproblemsK. WebbMAE 4020/5020

Solving ODE’s – General Aproach6 Have an ODE that is some function of the independent anddependent variables:, Numerical solutions amounts to approximatingStarting at some known initial condition, 0 , propagate thesolution forward in time:oris called the increment function Represents a slope, though not necessarily the slope at ,is the time step:K. WebbMAE 4020/5020

One‐Step vs. Multi‐Step Methods7 One‐step methods Useonly information at current value of(i.e., or ) to determine the increment function, , tobe used to propagate the solution forward to Collectively known as Runge‐Kutta methods We’ll focus on these exclusively in this course Multi‐step methods Useboth current and past values ofinformation about the trajectory of Improved accuracyK. Webbto provideMAE 4020/5020

8Euler’s MethodWe’ll first look at three specific Runge‐Kuttaalgorithms, before returning to a developmentof the Runge‐Kutta approach from a moregeneral perspective.K. WebbMAE 4020/5020

Euler’s Method9 Given an ODE of the form,approximate the solution,, using the formulawhere the increment function is the current derivative, That is, assume the slope of K. WebbUse the slope at,is constant forto extrapolate toMAE 4020/5020

Euler’s Method10 Euler’s methodformula: Increment function isthe current slope:K. WebbMAE 4020/5020

Euler’s Method ‐ Example11 Use Euler’s method tosolve.50.5given an initial conditionof03and a step size of0.5 True solution is:.K. Webb5 .MAE 4020/5020

Euler’s Method ‐ Example12K. WebbMAE 4020/5020

Euler’s Method ‐ Error13 Two types of truncation error:Local – error due to the approximation associated with thegiven method over a single time step Global – error propagated forward from previous time steps Total error is the sum of local and global error Representing the solution to the ODE as a Taylor seriesexpansion about, the solution atis:, ,2! ,!Where the remainder term is:K. WebbMAE 4020/5020

Euler’s Method ‐ Error14 Euler’s method is the Taylor series, truncated after thefirst derivative term For small enough , the error is dominated by the nextterm in the series, so Local error is proportional to Analysis of the global (i.e. propagated) error is beyondthe scope of this course, but the result is that globalerror is proportional toK. WebbMAE 4020/5020

Euler’s Method – Stability15 Euler’s method will result in error, but worse yet, it may be unstable Unstable if errors grow without boundConsider, for example, the following ODE:, The true solution decays exponentially to zero: Using Euler’s method, the solution is1 This solution will grow without bound if 1 K. Webb1, i.e. if2/If the step size is too large, solution blows upEuler’s method is conditionally stableMAE 4020/5020

Stability of Euler’s Method – Examples16K. WebbMAE 4020/5020

17K. WebbHeun’s MethodMAE 4020/5020

Heun’s Method18 Euler’s assumes a constant slope for the incrementfunction: Improve accuracy of the solution by using a moreaccurate slope estimate for Heun’s method first applies Euler’s method to predictthe value of at– the predictor equation: This value is then used to predict the slope atK. WebbMAE 4020/5020

Heun’s Method19 The increment function is the average of the slopeand the slope atat The next value ofequation:K. Webbis given by the correctorMAE 4020/5020

Heun’s Method – Summary20 Apply Euler’s – thepredictor equation – topredict Calculate slope at Compute average of thetwo slopes Use slope average topropagate the solutionforward to– thecorrector equationK. WebbMAE 4020/5020

Heun’s Method – Example21K. WebbMAE 4020/5020

Heun’s Method with Iteration22 Predictor equation:, Corrector equation:,,2 The corrector equation can be applied iteratively, providing arefined estimate of Iterate until approximate error falls below some stopping criterion 100%K. WebbMAE 4020/5020

Iterative Heun’s Method – Algorithm23, 1 While,, 100% 1Does not necessarily converge to the correct solution,though will converge to a finite valueK. WebbMAE 4020/5020

Iterative Heun’s Method – Example24K. WebbMAE 4020/5020

25K. WebbMidpoint MethodMAE 4020/5020

Midpoint Method26 The slope at themidpoint of a timeinterval used as theincrement function Provides a moreaccurate estimate ofthe slope across theentire time intervalK. WebbMAE 4020/5020

Midpoint Method27 Apply Euler’s method toapproximate at midpoint, 2Slope estimate at midpoint:, Midpoint slope estimate isincrement function,K. WebbMAE 4020/5020

Midpoint Method – Example28K. WebbMAE 4020/5020

One‐Step Methods – Error29MethodLocal ErrorGlobal ErrorEuler’sHeun’s (w/o iter.)MidpointK. WebbMAE 4020/5020

30K. WebbRunge‐Kutta MethodsMAE 4020/5020

Runga‐Kutta Methods31 Euler’s, Heun’s, and midpoint methods are specificcases of the broader category of one‐step methodsknown as Runge‐Kutta methodsRunge‐Kutta methods all have the same general form The increment function has the following form is the order of the Runge‐Kutta method K. WebbWe’ll see that Euler’s is a first‐order method, while Heun’sand midpoint are both second‐orderMAE 4020/5020

Runge‐Kutta Methods32 The increment function iswhere,,, ,, , The ’s, ’s, and ’s are constants Can see that Euler’s method is first‐order withK. WebbMAE 4020/5020

Runge‐Kutta Methods33 To determine values of ’s, ’s, and ’s: Setthe Runge‐Kutta formula equal to a Taylor series ofthe same order Equate coefficients An under‐determined system results Arbitrarily set one constant and solve for others Procedure is the same for all orders We’llstep through the derivation of the second‐orderRunge‐Kutta formulasK. WebbMAE 4020/5020

Second‐Order Runge‐Kutta Methods34 Second‐order Runge‐Kutta:(1)where,(2), (3)Second‐order Taylor series:,,!(4)where,K. Webb(5)MAE 4020/5020

Second‐Order Runge‐Kutta Methods35 Substituting (5) into (4), and recognizing thatthe Taylor series becomes, (6)!Next, represent (3) as a first‐order Taylor series It’s a function of two variables, for which the first‐orderTaylor series has the following formΔ , ,,,Δ,ΔΔ(7)Using (7), (3) becomes,K. Webb(8)MAE 4020/5020

Second‐Order Runge‐Kutta Methods36 Substituting (2) and (8) into (1),,, Now, set (9) equal to (6), the Taylor series,,, (9),,(10)Equating the coefficients in (10) gives three equations with fourunknowns:1(11)(12)(13)K. WebbMAE 4020/5020

Second‐Order Runge‐Kutta Methods37 We have three equations in four unknowns1(11)(12)(13) An under‐determined system Aninfinite number of solutions Arbitrarily set one constant –– to a certain valueand solve for the other three constants Different solution for each value of– a family ofsolutionsK. WebbMAE 4020/5020

– Heun’s Method38 Arbitrarily setand solve for the other constants, 1,,1The second‐order Runge‐Kutta formula becomes1212where,, ,This is Heun’s method,,2K. WebbMAE 4020/5020

– Midpoint Method39 Arbitrarily setand solve for the other constants0, 1,,The second‐order Runge‐Kutta formula becomeswhere,, 2,2This is the midpoint method,K. WebbMAE 4020/5020

Fourth‐Order Runge‐Kutta40 The most commonly used Runge‐Kutta method is the fourth‐order methodDerivation proceeds similar to that of the second‐order method Under‐determined system – family of solutionsMost common fourth‐order Runge‐Kutta method:1226where,1,2121,212, The increment function is a weighted average of four different slopesK. WebbMAE 4020/5020

4th‐Order Runge‐Kutta – Algorithm411.2.3.4.5.Calculate the slope at this is,Use to approximate/from . Calculate the slopehere this isUse to re‐approx./from . Calculate the slopehere this isUse to approx.from .Calculate the slope here thisisCalculate as a weightedaverage of the four slopesK. WebbMAE 4020/5020

Fourth‐Order Runge‐Kutta – Example42K. WebbMAE 4020/5020

43K. WebbSystems of EquationsMAE 4020/5020

Higher‐Order Differential Equations44 The ODE solution techniques we’ve looked at so farpertain to first‐order ODE’s Can be extended to higher‐order ODE’s by reducingto systems of first‐order equations An‐order ODE can be represented as a system offirst‐order ODE’s Solution method is applied to each equation at eachtime step before advancing to the next time step We’ll now revisit the fourth‐order quarter‐carexample from the first day of classK. WebbMAE 4020/5020

Fourth‐Order ODE – Example45 Recall the quarter‐car model from the introductorysection of this course Apply Newton’s second law to eachmass to derive the governing fourth‐order ODE Single 4th‐order equation, orTwo 2nd‐order equations0 Want to reduce to a system of fourfirst‐order ODE’s K. WebbPut into state‐space formMAE 4020/5020

Fourth‐Order ODE – Example460(1)(2) Reducing the ODE to a system of first‐order ODE’s is very similar torepresenting our system in state‐space form: The only difference being that we ultimately won’t actually representthe system in matrix formDefine a state vector of displacements and velocities:(3)K. WebbMAE 4020/5020

Fourth‐Order ODE – Example47 Rewrite (1) and (2) using the state variables definedin (3)0(4)(5) The state variable representation of the system is00K. Webb001001000 (6)MAE 4020/5020

Fourth‐Order ODE – Example48 Equation (6) clearly shows our system of four first‐orderODE’s Alternatively, could have derived the state‐space equationsdirectly (e.g. using a bond graph approach)In MATLAB, we’ll represent our system as ann‐dimensional function A vector of n functions:(7)(8)(9)(10)K. WebbMAE 4020/5020

Fourth‐Order ODE – Example49 In MATLAB, define theshown below ‐order system of ODE’s asAn ‐dimensional functionHere, the ODE function includes parameters (etc.) in addition to variables and, , Can create an anonymous function wrapper in the calling m‐file to allow for the passing of parametersK. WebbMAE 4020/5020

Fourth‐Order ODE – Example50 Basic formula remains the same Advance the solution to the next time step using the incrementfunctionNow, the output is the vector of states, and the incrementfunction is an ‐dimensional vectoror, ,,, ,,,,,, ,,,, ,Requires only a minor modification of the code written forfirst‐order ODE’s to accommodate ‐dimensional functionsK. WebbMAE 4020/5020

Fourth‐Order ODE – Example51 Often want to pass parameters (i.e. Input arguments in addition toand ) to the ODE function Two options (see Section 2: Programming with MATLAB notes): K. WebbInclude a varargin input argument in the ODE solver definitionUse an anonymous function wrapper for the ODE function, e.g.:MAE 4020/5020

Fourth‐Order ODE – Example52K. WebbMAE 4020/5020

Fourth‐Order ODE – Example53K. WebbMAE 4020/5020

Fourth‐Order ODE – Example54K. WebbMAE 4020/5020

55K. WebbSolving ODE’s in MATLABMAE 4020/5020

MATLAB’s ODE Solvers56 MATLAB has several ODE solvers Stiff ODE’s are those with a large range of eigenvalues – i.e. both very fastand very slow system poles ode45.m should usually be first choice for non‐stiff problemsNumerical solution is difficultFrom the MATLAB documentation:SolverStiffnessAccuracyWhen to useode45Non‐stiffMediumMost of the time. First choice.ode23Non‐stiffLowFor problems with crude error tolerances or for solvingmoderately stiff problems.ode113Non‐stiffLow to highFor problems with stringent error tolerances or forsolving computationally intensive problems.ode15sStiffLow to mediumIf ode45 is slow because the problem is stiff.ode23sStiffLowIf using crude error tolerances to solve stiff systems.K. WebbMAE 4020/5020

Solving ODE’s in MATLAB – ode45.m57[t,y] ode45(dydt,tspan,y0,options)dydt: handle to the ODE function – n‐dimensional tspan: vector of initial and final times – [ti,tf] y0: initial conditions – an n‐vector options: structure of options created with odeset.m t: column vector of time points y: solution matrix – length(t)n Syntax for all other solvers is identicalode45 uses an adaptive algorithm that uses fourth‐ andfifth‐order Runge‐Kutta formulas K. WebbVariable step sizeMAE 4020/5020

Fourth‐Order ODE – Example58K. WebbMAE 4020/5020

Passing Parameters as varargin59K. WebbMAE 4020/5020

Exercise – Solving ODE’s in MATLAB60 A simple pendulum of length is described by thefollowing second‐order ODEsinExercise This can be reduced to a system of two first‐orderODE’s:sin Define a function to describe this system of ODE’sWrite an m‐file that uses ode45.m to determineand plotandfor 010 K. Webb0.510 and 175 0Use odeset.m to set Reltol to different values(e.g. 10e-3 and 10e-6) and notice the effect onthe stability for175 MAE 4020/5020

3 Ordinary Differential Equations K. Webb MAE 4020/5020 Differential equations can be categorized as either ordinary or partialdifferential equations Ordinarydifferential equations (ODE's) - functions of a single independent variable Partial differential equations (PDE's) - functions of two or more independent variables

Related Documents: