New User Special Price Expires in

Let's log you in.

Sign in with Facebook


Don't have a StudySoup account? Create one here!


Create a StudySoup account

Be part of our community, it's free to join!

Sign up with Facebook


Create your account
By creating an account you agree to StudySoup's terms and conditions and privacy policy

Already have a StudySoup account? Login here

Intermediate Mechanics

by: April Ferry

Intermediate Mechanics PHYS 321

April Ferry

GPA 3.63


Almost Ready


These notes were just uploaded, and will be ready to view shortly.

Purchase these notes here, or revisit this page.

Either way, we'll remind you when they're ready :)

Preview These Notes for FREE

Get a free preview of these Notes, just enter your email below.

Unlock Preview
Unlock Preview

Preview these materials now for free

Why put in your email? Get access to more of this material and other relevant free materials for your school

View Preview

About this Document

Class Notes
25 ?




Popular in Course

Popular in Physics 2

This 53 page Class Notes was uploaded by April Ferry on Thursday October 15, 2015. The Class Notes belongs to PHYS 321 at New Mexico Institute of Mining and Technology taught by Staff in Fall. Since its upload, it has received 233 views. For similar materials see /class/223640/phys-321-new-mexico-institute-of-mining-and-technology in Physics 2 at New Mexico Institute of Mining and Technology.


Reviews for Intermediate Mechanics


Report this Material


What is Karma?


Karma is the currency of StudySoup.

You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 10/15/15
Chapter 7 Ordinary Differential Equa ons MATLAB has several different functions for the numerical solution of ordinary dif ferential equations This chapter describes the simplest of these functions and then compares all of the functions for efficiency accuracy and special features Stiffness is a subtle concept that plays an important role in these comparisons 71 Integrating Differential Equations The initial value problem for an ordinary differential equation involves nding a function yt that satisfies 59fmm together With the initial condition 9750 yo A numerical solution to this problem generates a sequence of values for the indepen dent variable to t1i i i and a corresponding sequence of values for the dependent variable y0y1 i i i so that each yn approximates the solution at tn yn mat n 01Hi Modern numerical methods automatically determine the step sizes hn tn1 in so that the estimated error in the numerical solution is controlled by a specified tolerance The Fundamental Theorem of Calculus gives us an important connection be tween differential equations and integrals zh MHWFWW fmmmw 1 2 Chapter 7 Ordinary Differential Equations We cannot use numerical quadrature directly to approximate the integral because we do not know the function ys and so cannot evaluate the integrand Nevertheless the basic idea is to choose a sequence of values of h so that this formula allows us to generate our numerical solution One special case to keep in mind is the situation where ft y is a function of t alone The numerical solution of such simple differential equations is then just a sequence of quadraturesi n1 y y fsds tn Through this chapter we frequently use dot notation for derivatives 7 dyt 7 d2yt y dt andy M 72 Systems of Equations Many mathematical models involve more than one unknown function and second and higher order derivatives These models can be handled by making yt a vector valued function of ti Each component is either one of the unknown functions or one of its derivatives The MATLAB vector notation is particularly convenient here or example the secondorder differential equation describing a simple har monic oscillator in 7m becomes two rstorder equations The vector yt has two components zt and its rst derivative Mt 9W W i i Using this vector the differential equation is 7 W i i 710 l 92W 91 t The MATLAB function de ning the differential equation has t and y as input arguments and should return ft y as a column vector For the harmonic oscillator the function is function ydot harmonic ty ydot y2 y1 A fancier version uses matrix multiplication in an inline function f inline 0 1 1 oy t y 73 Linearized Differential Equations 3 In both cases the variable t has to be included as the rst argument even though it is not explicitly involved in the differential equation slightly more complicated example the twobody problem describes the orbit of one body under the gravitational attraction of a much heavier bodyi Using Cartesian coordinates ut and vt centered in the heavy body the equations are 1W utTt3 at 7vlttgtrlttgt3 Where Tt xut2 vt2 The vector yt has four components W The differential equation is 1W 7 W 7ulttgtTlttgt3 vtTt3 The MATLAB function could be function ydot twobodyty r sqrty1 2 y2 2 ydot y3 y4 y1r 3 y2r 3 A more compact MATLAB function is function ydot twobodyty ydot y34 y12normy12 3 Despite the use of vector operations the second M le is not signi cantly more ef cient than the rst 73 Linearized Differential Equations The local behavior of the solution to a differential equation near any point tmyc can be analyzed by expanding ft y in a twodimensional Taylor series ay ftcyycatitcJyycm Where 5f 7 0 i at 0790 J 6y toyyc 4 Chapter 7 Ordinary Differential Equations The most important term in this series is usually the one involving J the Jacobiani For a system of differential equations With n components y1t f1tyyiy 7yn d 9205 7 f2lttvylvwwyn 3 z z ynt fnltt7y1v v 7 the Jacobian is an n byn matrix of partial derivatives i y1 y2 3y J 33 11 3112 3y 3111 3112 By The in uence of the Jacobian on the local behavior is determined by the solution to the linear system of ordinary differential equations y39Jy Let Ak Mk 21 be the eigenvalues of J and A diagk the diagonal eigenvalue matrixi If there is a linearly 39 A set of col quot 39 V then J VAV 1 The linear transformation VI y transforms the local system of equations into a set of decoupled equations for the individual components of x 13916 Akzk The solutions are mt eAktit ztc A single component 1k t grows With t if Mk is positive decays if Mk is negative and oscillates if lk is nonzeroi The components of the local solution yt are linear combinations of these behaviors For example the harmonic oscillator 701 9 7109 is a linear systemi The Jacobian is simply the matrix hi5 El 74 SingleStep Methods 5 The eigenvalues of J are ii and the solutions are purely oscillatory linear combi nations of e and 64 i A nonlinear example is the twobody problem 93 M t t y l mom y2tTt3 where W y1t2 y2t2 In an exercise we ask you to show that the Jacobian for this system is 0 0 T5 0 i 0 0 0 T5 r5 2y 7 yg 3y1y2 0 0 3y1y2 2y 7 y 0 0 J It turns out that the eigenvalues of J just depend on the radius Tt 1 239 T32 7239 A We see that one eigenvalue is real and positive so the corresponding component of the solution is growing One eigenvalue is real and negative corresponding to a decaying componenti Two eigenvalues are purely imaginary corresponding to os cillatory components However the overall global behavior of this nonlinear system is quite complicated and is not described by this local linearized analysis 74 SingleStep Methods The simplest numerical method for the solution of initial value problems is Euler s method It uses a xed step size h and generates the approximate solution by yn1 yn hftmyn tn1 in h The MATLAB code would use an initial point to a nal point tfinal an initial value yO a step size h and an inline function or function handle fl The primary loop would simply be while 1 lt tfinal y y hfeva1fty t t h end 6 Chapter 7 Ordinary Differential Equations Note that this works perfectly well if yO is a vector and f returns a vector As a quadrature rule for integrating ft Eulerls method corresponds to a rectangle rule where the integrand is evaluated only once at the left hand endpoint of the interval It is exact if is constant but not if is linear So the error is proportional to h Tiny steps are needed to get even a few digits of accuracy But from our point of view the biggest defect of Eulerls method is that it does not provide an error estimate There is no automatic way to determine what step size is needed to achieve a speci ed accuracyi lf Eulerls method is followed by a second function evaluation we begin to get a viable algorithmi There are two natural possibilities corresponding to the midpoint rule and the trapezoid rule for quadrature The mi point analog uses Euler to step halfway across the interval evaluates the function at this intermediate point then uses that slope to take the actual step 81 ftmyn h h 82 fan 5797 351 yn1 M h82 tn1 in h The trapezoid analog uses Euler to take a tentative step across the interval evaluates the function at this exploratory point then averages the two slopes to take the actual step 81 ftmyn s2 ftn h yn h51 8 3 yn1 yn h tn1 in h If we were to use both of these methods simultaneously they would produce two different values for yn1i The difference between the two values would provide an error estimate and a basis for picking the step size Furthermore an extrapolated combination of the two values would be more accurate than either one individuallyi Continuing with this approach is the idea behind singlestep methods for inte grating ODEsi The function ft y is evaluated several times for values oft between tn and tn1 and values of y obtained by adding linear combinations of the values of f to yni The actual step is taken using another linear combination of the function valuesi Modern versions of singlestep methods use yet another linear combination of function values to estimate error and determine step size Singlestep methods are often called RunyeKutta methods after the two Ger man applied mathematicians who rst wrote about them around 1905 The classical RungeKutta method was widely used for hand computation before the invention of digital computers and is still popular today It uses four function evaluations per step 81 ftmyn gt SingleStep Methods 7 h h 82 in Evyn l 551 h h 83 flttn 79n532 s4 flttnh7ynh53 h yn1 yngsl23223334 tn1 tnh If fty does not depend on y then classical RungeKutta has 32 33 and the method reduces to Simpson s quadrature ru er Classical RungeKutta does not provide an error estimate The method is sometimes used with a step size h and again with step size h2 to obtain an error estimate but we now know more ef cient methods Several of the ODE solvers in MATLAB including the textbook solver we describe later in this chapter are singlestep or Runge Kutta solversi A general singlestep method is characterized by a number of parameters ai BM 71 and 61x There are k stages Each stage computes a slope 31 by evaluating fty for a particular value oft and a value of y obtained by taking linear combinations of the previous slopesi 13971 3239 flttnaihvynhZ Lj3jgt7 i17k j1 The proposed step is also a linear combination of the slopes k yn1 M h Z Y i i1 An estimate of the error that would occur with this step is provided by yet another linear combination of the slopes k en1 h Z isi i1 If this error is less than the speci ed tolerance then the step is successful and yn1 is accepted f not the step is a failure and yn1 is rejected In either case the error estimate is used to compute the step size h for the next step The parameters in these methods are determined by matching terms in Taylor series expansions of the slopes These series involve powers of h and products of various partial derivatives of ft The order of a method is the exponent of the smallest power of h that cannot be matched It turns out that one two three or four stages yield methods of order one two three or four respectivelyi But it takes six stages to obtain a fthorder method The classical Runge Kutta method has four stages and is fourthorderi The names of the MATLAB ODE solvers are all of the form odennxx with digits 1111 indicating the order of the underlying method and a possibly empty xx indicating 8 Chapter 7 Ordinary Differential Equations some special characteristic of the method If the error estimate is obtained by comparing formulas with different orders the digits 1111 indicate these orders For example ode45 obtains its error estimate by comparing a fourthorder and a fth order formula 75 The 3523 Algorithm Our textbook function ode23tx is a simpli ed version of the function ode23 that is included with MATLAB The algorithm is due to Bogacki and Shampine 3 6 The 23 in the function names indicates that two simultaneous single step formulas one of second order and one of third order are involved method has three stages but there are four slopes 3139 because after the rst step the 31 for one step is the 34 from the previous step The essentials are 81 ftmyn h h 32 ftn 57yn 581 a 27 39D 3 3 33 ftn 1hwy 1hs2 tn tn h yn1 yn 3281 382 483 84 flttn1yyn1 n1 7227531 632 833 7 934 The simpli ed pictures in gure 71 show the starting situation and the three stages We start at a point tmyn with an initial slope 31 ftnyn and an estimate of a good step size h Our goal is to compute an approximate solution yn1 at tn1 tn h that agrees with true solution ytn1 to within the speci ed tolerancesi The rst stage uses the initial slope 31 to take an Euler step halfway across the interval The function is evaluated there to get the second slope 32 This slope is used to take an Euler step three quarters of the way across the interval The function is evaluated again to get the third slope 33 A weighted average of the three slopes 1 3 6231 332 433 is used for the nal step all the way across the interval to get a tentative value for yn1i The function is evaluated once more to get 34 The error estimate then uses all four slopesi h en1 57531 632 833 7 934 If the error is within the speci ed tolerance then the step is successful the tentative value of yn1 is accepted and 34 becomes the 31 of the next step If the error is U1 The 8523 Algorithm 9 tn tnh tn tnh2 s3 ynp1 tn tn3hl4 tn tnh Figure 71 3325 algorithm too large then the tentative yn1 is rejected and the step must be redone In either case the error estimate en1 provides the basis for determining the step size h for the next step The rst input argument of ode23tx speci es the function ft This argu ment can take any of three different forms 0 A character string involving 1 andor y 0 An inline function 0 A function handle The inline function or the M le should accept two arguments usually but not necessarily t and y The result of evaluating the character string or the function should be a column vector containing the values of the derivatives dydtt The second input argument of ode23tx is a vector tspan with two compo nents t0 and ti inalt The integration is carried out over the interval to S t S tfinal One of the simpli cations in our textbook code is this form of tspant Other MATLAB ODE solvers allow more exible speci cations of the integration interval The third input argument is a column vector yO providing the initial value of yo yt0t The length of y0 tells ode23tx the number of differential equations in the systemt 10 Chapter 7 Ordinary Differential Equations A fourth input argument is optional and can take two different forms The simplest and most common form is a scalar numerical value rtol to be used as the relative error tolerance The default value for rtol is 10 3 but you can provide a different value if you want more or less accuracy The more complicated possibility for this optional argument is the structure generated by the MATLAB function odeseti his function takes pairs of arguments that specify many different options for the MATLAB ODE solversi For ode23tx you can change the default values of three quantities the relative error tolerance the absolute error tolerance and the M le that is called after each successful step The statement opts odeset reltol 1e5 abstol 1e8 outputfcn meodeplot creates a structure that speci es the relative error tolerance to be 10 5 the absolute error tolerance to be 10 3 and the output function to be myodeploti The output produced by ode23tx can be either graphic or numeric With no output arguments the statement ode23txFtspany0 produces a dynamic plot of all the components of the solution With two output arguments the statement toutyout ode23txFtspany0 generates a table of values of the solution 76 ode23tx Letls examine the code for ode23txi Here is the preamble function toutyout ode23txFtspany0arg4varargin oAUDEZSTX Solve nonstiff differential equations quot0 Textbook version of ODE23 oo UDE23TXFTSPANY0 with TSPAN T0 TFINAL integrates quot0 the system of differential equations y fty from quot0 t To to t TFINAL The initial condition is quot0 yT0 YO The input argument F is the name of an quot0 Mfile or an inline function or simply a character quot0 string defining fty This function must have two quot0 input arguments t and y and must return a column vector quot0 of the derivatives y quot0 With two output arguments TY UDE23TX returns quot0 a column vector T and an array Y where Yk is the quot0 solution at Tk 76 ode23tx 11 X With no output arguments 0DE23TX plots the solution X ODE23TXFTSPANYORTUL uses the relative error Z tolerance RTUL instead of the default 1 e3 Z ODE23TXFTSPANYOUPIS where UPTS UDESET reltol Z RTUL abstol ATUL outputfcn PLTFN uses relative error Z RTUL instead of 1 e3 absolute error ATUL instead of Z 1 e6 and calls PLTFN instead of ODEPLOT after each step Z More than four input arguments ODE23TXFTSPANYORTUL Z P1P2 are passed on to F FTYP1P2 Z ODEZSTX uses the RungeKutta 23 method of Bogacki and Z Shampine Z Example Z tspan 0 2pi Z y0 1 0 Z F O 1 10y Z ode23txFtspany0 Z See also UDE23 Here is the code that parses the arguments and initializes the internal variables rtol atol plotfun odeplot if nargin gt 4 amp isnumericarg4 rtol arg4 dm nagn 4ampiUmHm if isemptyarg4RelTol rtol arg4 RelTol end if isemptyarg4AbsTol atol arg4 AbsTol end if isemptyarg40utputhn plotfun arg4 Uutputhn end end to tspan1 tfinal tspan2 tdir signtfinal t0 plotit nargout 0 threshold atol rtol hmax abs01tfinalt0 t to y y0 1e3 1e Make F callable by feval Chapter 7 Ordinary Differential Equations if ischarF amp existF 2 F inlineFrt y elseif isaF sym F inlinecharF t y end quot0 Initialize output if plotit feva1plotfuntspany init else tout t yout y end The computation of the initial step size is a delicate matter because it requires some knowledge of the overall scale of the problemi s1 feva1F t y varargin r norms1maxabsythresholdinf realmin h tdir08rtol 13r Here is the beginning of the main loop The integration starts at t to and increments t until it reaches tfinali It is possible to go backward that is7 have t final lt to while t tfinal hmin 16epsabst if absh gt hmax h tdirhmax end if absh lt hmin h tdirhmin end quot0 Stretch the step if t is close to tfinal if 11absh gt abstfina1 t h tfinal t end Here is the actual computation The rst slope s1 has already been computed The function de ning the differential equation is evaluated three more times to obtain three more slopesi s2 feva1F th2 yh2s1 varargin s3 feva1F t3h4 y3h4s2 varargin tnew t h ynew y h2s1 3s2 4sS9 s4 feva1F tnew ynew varargin 76 ode23tx 13 Here is the error estimate The norm of the error vector is scaled by the ratio of the absolute tolerance to the relative tolerance The use of the smallest floatingpoint number realmin prevents err from being exactly zeroi e h 5s1 6s2 8s3 9s472 norme maxmaxabsy absynew threshold inf realmin err Here is the test to see if the step is successful If it is the result is plotted or appended to the output vectori If it is not the result is simply forgotteni if err lt rtol t tnew Y YneW if plotit if fevalplotfunty break end else toutend11 t youtend1 y end s1 s4 quot0 Reuse final function value to start new step end The error estimate is used to compute a new step size The ratio rtolerr is greater than one if the current step is successful or less than one if the current step failsi A cube root is involved because the BS2 is a thirdorder method This means that changing tolerances by a factor of eight Will change the typical step size and hence the total number of steps by a factor of two The factors 08 and 5 prevent excessive changes in step size quot0 Compute a new step size h hmin508rtolerrquot13 Here is the only place Where a singularity would be detected if absh lt hmin warningsprintf Step size oe too small at t t tfinal end end That ends the main loop The plot function might need to nish its work oen ht if plotit feval plotfun l done end 14 Chapter 7 Ordinary Differential Equations 77 Examples Please sit down in front of a computer running MATLABi Make sure ode23tx is in your current directory or on your MATLAB path Start your session by entering ode23tx 0 0 10 1 This should produce a plot of the solution of the initial value problem dy i0 dt y01 09g 10 The solution of course is a constant function yt 1 Now you can press the up arrow key use the left arrow key to space over to the 0 and change it to something more interesting Here are some examples At rst welll change just the 0 and leave the O 10 and 1 alone F Exact solution 707 t 1th22 y expt y expt 11 3t 11og13t3 Singular 2y yquot2 21exp2t Make up some of your own examples Change the initial condition Change the accuracy by including 1e6 as the fourth argumenti Now let s try the harmonic oscillator a secondorder differential equation writ ten as a pair of two rstorder equationsi First create an inline function to specify the equations Use eit er F inline y2 y1 t y 01 11 inline 0 1 1 0y t y Then the statement ode23txF 0 2pi 1 0 plots two functions oft that you should recognize If you want to produce a phase plane plot you have two choices One possibility is to capture the output and plot it after the computation is complete ty ode23txF0 2pi1 0 ploty1y2 o axis12 12 12 121 axis square 77 Examples 15 The more interesting possibility is to use a function that plots the solution While it is being computedl MATLAB provides such a function in odephas2 m It is accessed by using odeset to create an options structure opts odeset re1tol 1e4 abstol 1e6 outputfcn odephas2 If you want to provide your own plotting function it should be something like function flag phaseplot tyjob persistent p if isequa1job init p ploty1 y2 o erasemode none axis12 12 12 121 axis square flag 0 elseif isequa1job set p xdata y1 ydata y2 This is With opts odeset re1tol 1e4 abstol 1e6 outputfcn phaseplot Once you have decided on a plotting function and created an options structure you can compute and simultaneously plot the solution With ode23txF 0 2pi 1 0 opts Try this With other values of the tolerancesl lssue the command type twobody to see if there is an M le twobody m on your path If not nd the three lines of code earlier in this chapter and create your own M lel Then try ode23tx twobody O 2pi 1 O 0 1 The code and the length of the initial condition indicate that the solution has four components But the plot shows only three Why Hint nd the zoom button on the gure Window toolbar and zoom in on the blue curvel ou can vary the initial condition of the twobody problem by changing the fourth componentl y0 1 0 O changethis ode23tx twobody O 2pi yO Graph the orbit and the heavy body at the origin With 16 Chapter 7 Ordinary Differential Equations y0 1 O 0 changethis ty ode23txtwobody 0 2pi yO ploty 1 y 2 00 ro axis equal You might also want to use something other than 27f for tf inali 78 Lorenz Attractor One of the world s most extensively studied ordinary differential equations is the Lorenz chaotic attractori It was first described in 1963 by Edward Lorenz an MIT mathematician and meteorologist who was interested in uid ow models of the earthls atmospherei An excellent reference is a book by Colin Sparrow We have chosen to express the Lorenz equations in a somewhat unusual way involving a matrixvector product 9 Ay The vector y has three components that are functions of t 91 Mt 92W 9W Despite the way we have written it this is not a linear system of differential equa tions Seven of the nine elements in the 3 by3 matrix A are constant but the other two depend on y2ti 5 0 92 A 0 70 a y2 P 1 The first component of the solution7 y1t is related to the convection in the atmo spheric ow while the other two components are related to horizontal and vertical temperature variation The parameter a is the Prandtl number p is the normal ized Rayleigh number and 6 depends on the geometry of the domain The most popular values of the parameters a 10 p 28 and 83 are outside the ranges associated with the earthls atmosphere The deceptively simple nonlinearity introduced by the presence of yg in the system matrix A changes everything There are no random aspects to these equa tions so the solutions yt are completely determined by the parameters and the initial conditions but their behavior is very difficult to predict For some values of the parameters the orbit of yt in threedimensional space is known as a strange attractori It is bounded but not periodic and not convergenti It never intersects itself It ranges chaotically back and forth around two different points or attractorsi For other values of the parameters the solution might converge to a xed point diverge to in nity or oscillate periodicallyi 78 Lorenz Attractor 17 Figure 73 Phase plane plot of Lorenz attractor Let7s think of 77 yg as a free parameter7 restrict p to be greater than one7 and study the matrix 6 0 77 A 0 70 a in P 1 It turns out that A is singular if and only if 77 iv w 1 18 Chapter 7 Ordinary Differential Equations The corresponding null vector normalized so that its second component is equal to 77 is p 7 1 n n With two different signs for 77 this de nes two points in threedimensional space These points are xed points for the differential equation If p 7 1 9750 77 n then for all t 0 W 0 0 and so yt never changes However these points are unstable xed points If yt does not start at one of these points it will never reach either of them if it tries to approach either point it will be repulsedi We have provided an M le lorenzgui m that facilitates experiments with the Lorenz equationsi Two of the parameters 6 83 and a 10 are xed A uicontrol offers a choice among several different values of the third parameter pi A simpli ed version of the program for p 28 would begin with rho 28 sigma 10 beta 83 eta sqrtbetarho1 A I beta 0 eta O sigma sigma eta rho 1 J The initial condition is taken to be near one of the attractorsi yc rho1 eta eta y0 yc 0 0 3 The time span is in nite so the integration will have to be stopped by another uicontroli tspan O Inf opts odeset re1tol 1e6 outputfcn 1orenzplot ode451orenzeqn tspan yO opts A The matrix A is passed as an extra parameter to the integrator which sends it on to lorenzeqn the subfunction de ning the differential equation The extra parameter machinery included in the function functions allows lorenzeqn to be written in a particularly compact manneri 79 Stiffness 19 function ydot lorenzeqntyA A13 y2 A31 y2 ydot Ay Most of the complexity of lorenzgui is contained in the plotting subfunction lorenzploti It not only manages the user interface controls it must also anticipate the possible range of the solution in order to provide appropriate axis scalingi 79 Stiffness Stiffness is a subtle dif cult and important concept in the numerical solution of ordinary differential equations It depends on the differential equation the initial conditions and the numerical method Dictionary de nitions of the word stiff involve terms like not easily bent rigid and stubborn We are concerned with a computational version of these properties A problem is sti the solution being sought is varying slowly but there are nearby solutions that vary rapidly so the numerical method must take small steps to obtain satisfactory results Stiffness is an ef ciency issue If we werenlt concerned with how much time a computation takes we wouldnlt be concerned about stiffness Nonstiff methods can solve stiff problems they just take a long time to do it A model of ame propagation provides an example We learned about this example from Larry Shampine one of the authors of the MATLAB ODE suitei If you light a match the ball of ame grows rapidly until it reaches a critical size Then it remains at that size because the amount of oxygen being consumed by the combustion in the interior of the ball balances the amount available through the surface The simple model is o 7y3 y05 0 t 26 The scalar variable yt represents the radius of the ball The y2 and y3 terms come from the surface area and the volume The critical parameter is the initial radius 6 which is small We seek the solution over a length of time that is inversely proportional to 6 At this point we suggest that you start up MATLAB and actually run our examples It is worthwhile to see them in action We will start with ode45 the workhorse of the MATLAB ODE suitei If 6 is not very small the problem is not very stiff Try 6 01 and request a relative error of 10 4i delta 001 F inline y 2 y 3 t y opts odeset Re1Tol 1e4 ode45F 0 2de1ta deltaopts 20 Chapter 7 Ordinary Differential Equations With no output arguments ode45 automatically plots the solution as it is computed You should get a plot of a solution that starts at y 01 grows at a modestly increasing rate until t approaches 100 which is 16 then grows rapidly until it reaches a value close to l where it remains Now let s see stiffness in action Decrease 6 by a couple of orders of magnitude If you run only one example run this one delta 00001 ode45F 0 2de1ta deltaopts ode45 i 100017 Figure 74 375239 behavior of ode45 You should see something like gure 74 although it will take a long time to complete the plot If you get tired of watching the agonizing progress click the stop button in the lower left corner of the windowi Turn on zoom and use the mouse to explore the solution near where it rst approaches steady state You should see something like the detail in the gure 74 Notice that ode45 is doing its job ltls keeping the solution within 10 4 of its nearly constant steady state value But it certainly has to work hard to do it If you want an even more dramatic demonstration of stiffness decrease the tolerance to 1075 or 1075 This problem is not stiff initiallyi It only becomes stiff as the solution ap proaches steady state This is because the steady state solution is so rigidi Any solution near yt 1 increases or decreases rapidly toward that solution We 79 Stiffness 21 should point out that rapidly here is with respect to an unusually long time scale What can be done about stiff problems You donlt want to change the dif ferential equation or the initial conditions so you have to change the numerical method Methods intended to solve stiff problems ef ciently do more work per step but can take much bigger steps Stiff methods are implicit At each step they use MATLAB matrix operations to solve a system of simultaneous linear equations that helps predict the evolution of the solution For our ame example the matrix is only lbyl but even here stiff methods do more work per step than nonstiff methods odeZas 10001 7 l l l l 098 1 102 104 106 108 11 112 Figure 75 Sti behavior of ode235 Letls compute the solution to our ame example again this time with one of the ODE solvers in MATLAB whose name ends in s for stiff delta 00001 ode235 F O 2de1tade1taopts Figure 75 shows the computed solution and the zoom detail You can see that ode235 takes many fewer steps than ode45 This is actually an easy problem for a stiff solver In fact ode235 takes only 99 steps and uses just 412 function evalua tions while ode45 takes 3040 steps and uses 20179 function evaluations Stiffness even affects graphical output The print files for the ode45 gures are much larger than those for the ode235 gures 22 Chapter 7 Ordinary Differential Equations Imagine you are returning from a hike in the mountains You are in a narrow canyon with steep slopes on either side An explicit algorithm would sample the local gradient to nd the descent direction But following the gradient on either side of the trail will send you bouncing back and forth across the canyon as with ode45i You will eventually get home but it will be long after dark before you arrive An implicit algorithm would have you keep your eyes on the trail and anticipate where each step is taking you It is well worth the extra concentration This ame problem is also interesting because it involves something called the Lambert W function The differential equation is separablei Integrating once gives an implicit equation for y as a function of ti 1 1 1 1 Zlogltgilgt 3loglt371gt7t This equation can be solved for y The exact analytical solution to the ame model turns out to be 1 W where a 16 7 1 The function Wz the Lambert W function is the solution to WzeWz 2 With MATLAB and the Symbolic Math Toolbox connection to Maple the statements y dsolve Dy y 2 y 3 y0 1100 simplifyy pretty y ezplot y 0 200 amplt produce 1ambertw99 exp99 11 1 and the plot of the exact solution shown in gure 716 If the initial value 1100 is decreased and the time span 0 S t S 200 increased the transition region becomes narroweri The Lambert W function is named after Jr Hi Lambert 1728 7 1777 Lambert was a colleague of Euler and Lagrange at the Berlin Academy of Sciences and is best known for his laws of illumination and his proof that 7r is irrational The function was rediscovered a few years ago by Corless Gonnet Hare and Jeffrey working on Maple and by Don Knuth 710 Events So far we have been assuming that the tspan interval to S t S tfinal is a given part of the problem speci cation or we have used an in nite interval and a GUI 710 Events 23 l l l l 0 20 40 60 80 100 120 140 160 180 200 t Figure 76 Exact solution for the flame example button to terminate the computation In many situations the determination of tfinal is an important aspect of the probleml example is a body falling under the force of gravity and encountering air resistance When does it hit the ground Another example is the twobody problem the orbit of one body under the gravitational attraction of a much heavier body What is the period of the orbit The events feature of the MATLAB ODE solvers provides answers to such questions Events detection in ordinary differential equations involves two functions ft y and gt y and an initial condition to yo The problem is to nd a function yt and a nal value t so that y WM 9050 yo and yOnWQ 0 A simple model for the falling body is y 71 y with initial conditions y0 l 0 The question is for what t does yt 0 The code for the function ft y is function ydot fty ydot y2 1y2 2 With the differential equation written as a rstorder system y becomes a vector with two components and so gt y yll The code for gty is function gstopisterminaldirection gty gstop y1 isterminal 1 direction 1 24 Chapter 7 Ordinary Differential Equations The rst output gstop is the value that we want to make zeroi Setting the second output isterminal to one indicates that the ODE solver should terminate when gstop is zero Setting the third output direction to the empty matrix indicates that the zero can be approached from either direction With these two functions available the following statements compute and plot the trajectory opts odeset events g y0 1 0 tytfina1 ode45f0 Infy0opts tfinal plotty1 O tfinal 1 0 o axis1 tfina11 1 1 1 x1abe1 t ylabe1 y title Fa11ing body text12 O tfinal num2strtfina1 The terminating value of t is found to be tfinal 16585 Falling body t nal16585 Figure 77 Event handling for falling object The three sections of code for this example can be saved in three separate M les with two functions and one script or they can all be saved in one function M lei In the latter case f and g become subfunctions and have to appear after the main body of code Events detection is particularly useful in problems involving periodic phenom ena The twobody problem provides a good example Here is the rst portion of a function M le orbit mi The input parameter is reltol the desired local relative tolerance function orbit relt 01 710 Events y0 1 O 0 03 opts odeset events Qgstop reltol reltol tyteye ode45twobody0 2piy0optsy0 tfinal teend yfinal yeend12 ploty1y2 00 ro axis1 105 35 351 The function ode45 is used to compute the orbit The rst input argument is a function handle twobody that references the function de ning the differential equations The second argument to ode45 is any overestimate of the time interval required to complete one period The third input argument is yO a 4vector that provides the initial position and velocity The light body starts at 10 which is a point with a distance 1 from the heavy body and has initial velocity 003 which is perpendicular to the initial position vector The fourth input argument is an options structure created by odeset that overrides the default value for relt 01 and that speci es a function gstop that de nes the events we want to locate The last argument is yO an extra argument that ode45 passes on to both twobody and gstopi The code for twobody has to be modi ed to accept a third argument even though it is not used function ydot twobodytyy0 r sqrty1 2 y2 2 ydot y3 y4 y1r 3 y2r 3 The ODE solver calls the gstop function at every step during the integration This function tells the solver whether or not it is time to stop function valistermdir gstoptyy0 d y12y012 v y34 val d v isterm 1 dir 1 The 2vector d is the difference between the current position and the starting point The 2vector v is the velocity at the current position The quantity val is the inner product between these two vectors Mathematically the stopping function is 9t7y WWW where d 9105 9107y2t 920T Points where gt 0 are the local minimum or maximum of dti By setting dir 1 we indicate that the zeros of gty must be approached from above so they correspond to minimal By setting isterm 1 we indicate that computation 26 Chapter 7 Ordinary Differential Equations of the solution should be terminated at the rst rninirnurni If the orbit is truly periodic7 then any minima of d occur When the body returns to its starting point al ing orbit With a very loose tolerance orbit20e3 produces tf inal 2 35087197761898 yf inal 0 98107659901079 0 00012519138559 and plots n o a i a 2 a 3 ti n 2 EI A nits Elia i Figure 78 Pen39odz39c orbit computed with loose tolerance You can see from both the value of yf inal and the graph that the orbit does not quite return to the starting point We need to request more accuracyi orbit10e6 produces tf inal 2 38025846171805 yf inal 0 99998593905521 0 00000000032240 Now the value of yfinal is close enough to yO that the graph of the orbit is effec tively closedi 711 Multistep Methods 27 711 Multistep Methods A singlestep numerical method has a short memory The only information passed from one step to the next is an estimate of the proper step size and perhaps the value of ftn yn at the point the two steps have in common As the name implies a multistep method has a longer memoryi After an initial startup phase a pth order multistep method saves up to perhaps a dozen values of the solution ynp1 ynp2i i yn71 yn and uses them all to compute yn1i In fact these methods can vary both the order p and the step size h Multistep methods tend to be more ef cient than singlestep methods for problems with smooth solutions and high accuracy requirements For example the orbits of planets and deep space probes are computed with multistep methodsi 712 The MATLAB ODE Solvers This section is derived from the Algorithms portion of the MATLAB Reference Man ual page for the ODE solversi ode45 is based on an explicit Runge Kutta 45 formula the Dormand Prince pairi It is a onestep solveri ln computing ytn1 it needs only the solution at the immediately preceding time point In general ode45 is the rst function to try for most problemsi ode23 is an implementation of an explicit Runge Kutta 23 pair of Bogacki and Shampinei It is often more ef cient than ode45 at crude tolerances and in the presence of moderate stiffnessi Like ode45 ode23 is a onestep solveri ode113 is a variableorder AdamsBashforth Moulton PECE solveri It is often more ef cient than ode45 at stringent tolerances and if the ODE le function is particularly expensive to evaluate ode113 is a multistep solver 7 it normally needs the solutions at several preceding time points to compute the current solution The above algorithms are intended to solve nonstiff systems If they appear to be unduly slow try using one of the stiff solvers belowi odelEs is a variableorder solver based on the numerical differentiation formu las NDFs Optionally it uses the backward differentiation formulas BDFs also known as Gearls method that are usually less ef cient Like ode113 odelEs is a multistep solveri Try odelEs if ode45 fails or is very inef cient and you suspect that the problem is stiff or if solving a differentialalgebraic problemi ode235 is based on a modi ed Rosenbrock formula of order 2 Because it is a onestep solver it is often more ef cient than odelEs at crude tolerancesi It can solve some kinds of stiff problems for which odelEs is not effective ode23t is an implementation of the trapezoidal rule using a free interpolanti Use this solver if the problem is only moderately stiff and you need a solution without numerical dampingi ode23t can solve DAEsi ode23tb is an implementation of TR BDFZ an implicit Runge Kutta formula with a rst stage that is a trapezoidal rule step and a second stage that is a backward differentiation formula of order two By construction the same iteration matrix is used in evaluating both stages Like ode235 this solver is often more ef cient than odelEs at crude tolerancesi 28 Chapter 7 Ordinary Differential Equations ere is a summary table from the MATLAB Reference Manual For each function it lists the appropriate problem type the typical accuracy of the method and the recommended area of usage 0 ode45i Nonstiff problems medium accuracyi Use most of the time This should be the rst solver you try 0 ode23i Nonstiff problems low accuracyi Use for large error tolerances or moderately stiff problems 0 ode113i Nonstiff problems low to high accuracy Use for stringent error tolerances or computationally intensive ODE functions 0 odelEsi Stiff problems low to medium accuracyi Use if ode45 is slow stiff systems or there is a mass matrix 0 ode235i Stiff problems low accuracy Use for large error tolerances with stiff systems or with a constant mass matrix 0 ode23ti Moderately stiff problems low accuracy Use for moderately stiff problems where you need a solution without numerical dampingi o ode23tbi Stiff problems low accuracyi Use for large error tolerances with stiff systems or if there is a mass matrix 713 Errors Errors enter the numerical solution of the initial value problem from two sources 0 Discretization error 0 Roundoff error Discretization error is a property of the differential equation and the numerical method If all the arithmetic could be performed with in nite precision discretiza tion error would be the only error presenti Roundoff error is a property of the computer hardware and the programi It is usually far less important than the discretization error except when we try to achieve very high accuracy Discretization error can be assessed from two points of view local and global Local discretization error is the error that would be made in one step if the previous values were exact and if there were no roundoff error Let unt be the solution of the differential equation determined not by the original initial condition at to but by the value of the computed solution at tni That is unt is the function oft de ned by an ft7un un n yn 713 Errors 29 The local discretization error dn is the difference between this theoretical solution and the computed solution ignoring roundoff determined by the same data at tni dn yn1 untn1 Global discretizatz39on error is the difference between the computed solution still ignoring roundoff and the true solution determined by the original initial con dition at to that is en yn 7 ytn The distinction between local and global discretization error can be easily seen in the special case where fty does not depend on y In this case the solution is simply an integral yt f739d739i Eulerls method becomes a scheme for numerical quadrature that might be called the composite lazy man7s rectangle rule77 It uses function values at the lefthand ends of the subintervals rather than at the midpoints EN N71 t frdr m E hum The local discretization error is the error in one subinterval n1 dn MW e frgtdr tn and the global discretization error is the total error N71 EN eN 2 MW e fTdT n0 0 In this special case each of the subintegrals is independent of the others the sum could be evaluated in any order so the global error is the sum of the local errors N71 EN Z dn 710 In the case of a genuine differential equation where fty depends on y the error in any one interval depends on the solutions computed for earlier intervalsi Consequently the relationship between the global error and the local errors is related to the stability of the differential equation For a single scalar equation if the partial derivative BfBy is positive then the solution yt grows as t increases and the global error will be greater than the sum of the local errors If BfBy is negative then the global error will be less than the sum of the local errors If BfBy changes sign or if we have a nonlinear system of equations where BfBy is a varying matrix the relationship between EN and the sum of the dn can be quite complicated and unpredictable Think of the local discretization error as the deposits made to a bank account and the global error as the overall balance in the account The partial derivative 30 Chapter 7 Ordinary Differential Equations BfBy acts like an interest rate If it is positive the overall balance is greater than the sum of the deposits If it is negative the final error balance might well be less than the sum of the errors deposited at each step Our code ode23tx like all the production codes in MATLAB only attempts to control the local discretization errori Solvers that try to control estimates of the global discretization error are much more complicated are expensive to run and are not very successful ndamental concept in assessing the accuracy of a numerical method is its orderi The order is defined in terms of the local discretization error obtained if the method is applied to problems with smooth solutions A method is said to be of order p if there is a number C so that idni S Chffl The number C might depend on the partial derivatives of the function defining the differential equation and on the length of the interval over which the solution is sought but it should be independent of the step number n and the step size hni The above inequality can be abbreviated using bigoh notation dn Oh 1 For example consider Eulerls method yn1 yn hnftm yn Assume the local solution un t has a continuous second derivative Then using Taylor series near the point tn un WW t tnuztn 0t 7inf Using the differential equation and the initial condition defining un t untn1 yn hnftmyn 0hi Consequently dn yn1 untn1 0013 We conclude that p 1 so Eulerls method is first order The MATLAB naming conventions for ODE solvers would imply that a function using Eulerls method by itself with fixed step size and no error estimate should be called Odell Now consider the global discretization error at a fixed point t t i As ac curacy requirements are increased the step sizes hn will decrease and the total number of steps N required to reach tf will increase Roughly we shall have if 7 to N 7 h where h is the average step size Moreover the global error EN can be expressed as a sum of N local errors coupled by factors describing the stability of the equations These factors do not depend in a strong way on the step sizes and so we can say 713 Errors 31 roughly that if the local error is Ohp1 then the global error will be NOhp1 0hpi This is why 10 1 was used instead of p as the exponent in the de nition of order For Eulerls method p 1 so decreasing the average step size by a factor of 2 will decrease the average local error by a factor of roughly 2p1 4 but about twice as many steps will be required to reach tf so the global error will be decreased by a factor of only 21 2 With higher order methods the global error for smooth solutions would be reduced by a much larger factor It should be pointed out that in discussing numerical methods for ordinary differential equations the word order can have any of several different meanings The order of a differential equation is the index of the highest derivative appearing For example dedtQ 7y is a secondorder differential equation The order of a system of equations sometimes refers to the number of equations in the systemi For example y39 2y 7 y22 72 yz is a secondorder systemi The order of a numerical method is what we have been discussing here It is the power of the step size that appears in the expression for the global errori One way of checking the order of a numerical method is to examine its behavior if fty is a polynomial in t and does not depend on y If the method is exact for tpil but not for tp then its order is not more than p The order could be less than p if the methods behavior for general functions does not match its behavior for polynomials Eulerls method is exact if ft y is constant but not if ft y t so its order is not greater than one With modern computers using IEEE oatingpoint doubleprecision arith metic the roundoff error in the computed solution only begins to become important if very high accuracies are requested or the integration is carried out over a long intervali Suppose we integrate over an interval of length L tf 7 to If the roundoff error in one step is of size 6 then the worst the roundoff error could be after N steps of size h is something like Le N e i h For a method with global discretization error of size Chp the total error would be something like Le 17 Chh For the roundoff error to be comparable with the discretization error we would nee 1 Le W h N lt6 The number of steps taken with this step size is roughly 1 C 71 N N LltEgt Here are the numbers of steps for various orders p if L 20 C 100 and e 2 32 Chapter 7 Ordinary Differential Equations p N 1 4 5 1017 3 5647721 5 37285 10 864 These values ofp are the orders for Eulerls method7 for the MATLAB functions ode23 and ode457 and a typical choice for the order in the variableorder method used by ode113i We see that the loworder methods have to take an impracti cally large number of steps before this worstcase roundoff error estimate becomes signi canti Even more steps are required if we assume the roundoff error at each step varies randomlyi The variableorder multistep function ode113 is capable of achieving such high accuracy that roundoff error can be a bit more signi cant with it 714 Performance We have carried out an experiment to see how all this applies in practice The differential equation is the harmonic oscillator 9W IO with initial conditions 10 lz390 07 over the interval 0 S t S 10 The interval is ve periods of the periodic solution7 so the global error can be computed simply as the difference between the initial and nal values of the solution Since the solution neither grows nor decays with t7 the global error should be roughly proportional to the local error The following MATLAB script uses odeset to change both the relative and the absolute tolerancesi The re nement level is set so that one step of the algorithm generates one row of output tol 10 k opts odeset re1tol tol abstol tol refine 1 tic ty ode23harmonic 0 10pi yO opts time toc steps 1engtht1 err maxabsyend y0 end The differential equation is de ned in harmonic mi function ydot harmonicty ydot y2 y1 714 Performance 33 error time Figure 79 Performance of ode solvers The script was run three times7 With ode237 ode457 and ode113l The rst plot in gure 79 shows how the global error varies With the requested tolerance for the three routinesl We see that the actual error tracks the requested tolerance quite well For ode23 the global error is about 36 times the tolerance7 for ode45 it is about 4 times the tolerance and for ode113 it varies between one and 45 times 34 Chapter 7 Ordinary Differential Equations the tolerance The second plot in gure 79 shows the numbers of steps required The results also t our model quite well Let 739 denote the tolerance i For ode23 the number of steps is about 107 13 Which is the expected behavior for a thirdorder method For ode45 the number of steps is about 97 15 Which is the expected behavior for a fthorder method For ode113 the number of steps re ects the fact that the solution is very smooth so the method was often able to use its maximum order 13 The third plot in gure 79 shows the execution times in seconds on a 800 MHZ Pentium III laptopi For this problem ode45 is the fastest method for tol erances of roughly 10 6 or larger While ode113 is the fastest method for more stringent tolerancesi The loworder method ode23 takes a very long time to ob tain high accuracy This is just one experiment on a problem With a very smooth and stable solution 715 Further Reading The MATLAB ODE suite is described in Additional material on the numerical solution of ordinary differential equations and especially stiffness is available in Ascher and Petzold l Brennan Campbell and Petzold 2 and Shampine Exercises 7ili a Show experimentally or algebraically that the BS2 algorithm is exact for my 1 my t and my K but not for my t b When is the ode23 error estimator exact 72 The error function erfz is usually de ned by an integral 1 gt 2 f er I 7 e I V77 0 but it can also be de ned as the solution to the differential equation 41 y0 0 Use ode23tx to solve this differential equation on the interval 0 lt z lt 2 Compare the results With the builtin MATLAB function erf x at the poiints chosen by ode23txi 73 a Write an M le named myrk4m in the style of ode23txm that imple ments the classical RungeKutta xed step size algorithmi Instead of an optional fourth argument rtol or opts the required fourth argument should be the step size h Here is the proposed preamble Exercises 7 5 X function toutyout myrk4FtspanyOhvarargin Z MYRK4 Classical fourthorder RungeKutta Z Usage is the same as 0DE23TX except the fourth argument Z is a fixed step size h K MYRK4FTSPANYOH with TSPAN TO IT integrates the Z system of differential equations y fty from t TO K to t IT The initial condition is yTO YO X With no output arguments MYRK4 plots the solution X With two output arguments TY MYRK4 returns T Z and Y so that Y k is the approximate solution at Tk Z More than four input arguments MYRK4 P1P2 Z are passed on to F FTYP1P2 b Roughly how should the error behave if the step size h for classical Runge Kutta is cut in half Hint Why is there a 477 in the name of myrk4 Run an experiment to illustrate this behavior c If you integrate the simple harmonic oscillator 7y over one full period 0 S t S 2 you can compare the initial and nal values of y to get a measure of the global accuracy If you use your myrk4 with a step size h 7r50 you should nd that it takes 100 steps and computes a result with an error about 10 Compare this with the number of steps required by ode23 ode45 and ode113 if the relative tolerance is set to 10 6 and the re nement level is set to one This is a problem with a very smooth solution so you should nd that ode23 requires more steps while ode45 and ode113 require fewer The ODE problem y39 71000y 7 sint cost y0 1 on the interval 0 S t S l is mildly stiff a Find the exact solution either by hand or using dsolve from the Symbolic Toolboxi b Compute the solution with ode23txi How many steps are required c Compute the solution with the stiff solver ode23si How many steps are required d Plot the two computed solutions on the same graph with line style for the ode23tx solution and o for the ode23s solutioni e Zoom in or change the axis settings to show a portion of the graph where the solution is varying rapidly You should see that both solvers are taking small steps f Show a portion of the graph where the solution is varying slowly You should see that ode23tx is taking much smaller steps than ode23si The following problems all have the same solution on 0 S t S 7r2i y39 cost y0 0 y39 1927 900 7y 90 07 90 1 36 Chapter 7 Ordinary Differential Equations 9 Sinty 90 07 90 1 a What is the common solution yt b Two of the problems involve second derivatives Rewrite these problems as rstorder systems y39 ft y involving vectors y and c What is the Jacobian J y for each problem What happens to each Jacobian as t approaches 7r2 d The work required by a Runge Kutta method to solve an initial value problem y39 fty depends on the function fty not just the solution Use odeset to set both reltol and abstol to 10 How much work does ode45 require to solve each problem Why are some problems more work than others e What happens to the computed solutions if the interval is changed to 0 S t S 7r f What happens on 0 S t S 7r if the second problem is changed to y39 i192i7 y0 0 76 Use the jacobian and eig functions in the Symbolic Toolbox to verify that the Jacobian for the twobody problem is 0 0 7 5 0 J 7 i 0 0 0 7 5 r5 2y 7 y 3y1y2 0 0 3y1y2 2y 7 y 0 0 and that its eigenvalues are 1 239 A T32 7i 77 a Verify that the matrix in the Lorenz equations 6 0 77 A 0 70 a in P 1 is singular if and only if 77 i 503 i 1 Verify that the corresponding null vector is p 7 l 77 77 b The Jacobian matrix J for the Lorenz equations is not A but is closely related to A Find J compute its eigenvalues at one of the xed points and verify that the xed point is unstable Exercises 37 78 5 to What happens if the initial value in lorenzgui is chosen to be very close to a critical point Roughly how far away from the critical point does the initial value have to be so that the trajectory starts towards the other critical point in a reasonable amount of time Does the answer depend on the value of p i All the values of p available with lorenzgui except p 28 give trajectories that eventually settle down to stable periodic orbits In his book on the Lorenz equations Sparrow classifies a periodic orbit by what we might call its signature a sequence of 7s and 77s specifying the order of the critical points that the trajectory circles during one period A single or 7 woul be the signature of a trajectory that circles just one critical point except that no such orbits exist The signature 77 indicates that the trajectory circles each critical point once The signature 7 7 777 would indicate a very fancy orbit that circles the critical points a total of eight times before repeating itself a What are the signatures of the four different periodic orbits generated by lorenzgui Be careful 7 each of the signatures is different and p 9965 is particularly delicatei b Modify lorenzgui to provide other values of pi Find a p that gives a periodic orbit with a different signature from the ones already generatedi c Find a p besides p 28 that gives a chaotic trajectoryi What are the periods of the periodic orbits generated for the different values of p available with lorenzgui The MATLAB demos directory contains an Mfile orbitode that uses ode45 to solve an instance of the restricted threebody problemi This involves the orbit of a light object around two heavier objects such as an Apollo capsule around the earth and the moon Run the demo and then locate its source code with the statements orbitode which orbitode Make your own copy of orbitodemi Find these two statements tspan O 7 y0 12 o o 1o4935750983031990726 These statements set the time interval for the integration and the initial position and velocity of the light objecti Our question is where do these values come from To answer this question find the statement 1 yte ye ie ode45 f tspany0 options Remove the semicolon and insert three more statements after it te ye ie 38 Chapter 7 Ordinary Differential Equations Run the demo againi Explain how the values of te ye and ie are related to tspan and in A classical model in mathematical ecology is the Lotka Volterra predator prey modeli Consider a simple ecosystem consisting of rabbits that have an in nite supply of food and foxes that prey on the rabbits for their food This is modeled by a pair of nonlinear firstorder differential equations d7 27 7an T0 7 0 d di 7fan7 f0 for where t is time Tt is the number of rabbits is the number of foxes and a is a positive constant If a 0 the two populations do not interact the rabbits do what rabbits do best and the foxes die off from starvation If a gt 0 the foxes encounter the rabbits with a probability that is proportional to the product of their numbers Such an encounter results in a decrease in the number of rabbits and for less obvious reasons an increase in the number of foxesi The solutions to this nonlinear system cannot be expressed in terms of other known functions the equations must be solved numerically It turns out that the solutions are always periodic with a period that depends on the initial conditions In other words for any T0 and f0 there is a value t tp when both populations return to their original valuesi Consequently for all 757 Tt tp TO ft tp W a Compute the solution with 7 0 300 f0 150 and a 001 You should nd that tp is close to 5 Make two plots one of T and f as functions oft and one a phase plane plot with T as one axis and f as the other b Compute and plot the solution with 7 0 15 f0 22 and a 001 You should find that tp is close to 662 c Compute and plot the solution with 7 0 102 f0 198 and a 001i Determine the period tp either by trial and error or with an event handler d The point To f0 la 2a is a stable equilibrium point If the popu lations have these initial values they do not change If the initial populations are close to these values they do not change very much Let ut 7 t71Ot and vt 7 2ai The functions ut and vt satisfy another nonlinear system of differential equations but if the uv terms are ignored the system becomes lineari What is this linear system What is the period of its periodic solutions Many modi cations of the Lotka Volterra predatorprey model see previous problem have been proposed to more accurately reflect what happens in nature For example the number of rabbits can be prevented from growing inde nitely by changing the first equation as follows dT T a 21 7 an T0T0 Exercises df if we ylt0gt 90 dt Where t is time Tt is the number of rabbits is the number of foxes a is a positive constant and R is a positive constant Because a is positive is negative Whenever 7 2 R Consequently the number of rabbits can never exceed R For a 0 01 compare the behavior of the original model With the behavior of this modi ed model With 400 In making this comparison solve t e equations With 7 0 300 and f0 150 over 50 units of time Make four different plots 0 Number of foxes and number of rabbits vs time for the original model 0 Number of foxes and number of rabbits vs time for the modi ed model 0 Number of foxes vsi number of rabbits for the original model 0 Number of foxes vsi number of rabbits for the modi ed model For all plots label all curves and all axes and put a title on the plot For the last two plots set the aspect ratio so that equal increments on the x and y axes are equal in size i An 80 kilogram paratrooper is dropped from an airplane at a height of 600 meters After 5 seconds the chute opens function of time yt is given by The paratrooper7s height as a Q 79atm y0 600 ft 90 0fts Where g 981 ms2 is the acceleration due to gravity and m 80 kg is the paratrooper7s mass The air resistance at is proportional to the square of the velocity With different proportionality constants before and after the chute opens at Kly39t2 tlt 5 seconds Kgy39t2 t 2 5 seconds a Find the analytical solution for the case K1 0 K2 0 At What height does the chute open How long does it take to reach the ground What is the impact velocity Plot the height vs time and label the plot appropriately b Consider the case K1 1150K2 4150 At What height does the chute open How long does it take to reach the ground What is the impact velocity Make a plot of the height vs time and label the plot appropriately i Determine the trajectory of a spherical cannonball in a stationary Cartesian coordinate system that has a horizontal zaxis a vertical yaxis and an origin at the launch point The initial velocity of the projectile in this coordinate system has magnitude v0 and makes an angle With respect to the zaxis 40 Chapter 7 Ordinary Differential Equations H of 90 radians The only forces acting on the projectile are gravity and the aerodynamic drag D which depends on the projectilels speed relative to any wind that might be present The equations describing the motion of the projectile are z39vcost9 y39vsini9 igcosi i v 79 7 g sini9 m Constants for this problem are the acceleration of gravity 9 9 81ms2 the mass m 15 kg and the initial speed v0 50 ms The wind is assumed to be horizontal and its speed is a speci ed function of time The aerodynamic drag is proportional to the square of the projectilels velocity relative to the wind Dt cps 7 1 10052 92 where c 02 is the drag coef cient p 1129 kgm3 is the density of air and s 0 25m2 is the projectilels cross sectional area Consider four different wind conditions 0 No windi wt 0 for all t o Steady head windi wt 710 ms for all t o lntermittent tail windi wt 10 ms if the integer part oft is even zero otherwise 0 Gusty windi wt is a Gaussian random variable with mean zero and standard deviation 10 ms The integer part of a real number t is denoted by it and is computed in MAT LAB by floor it A Gaussian random variable with mean 0 and standard a is generated by sigmarandn see the chapter on random numbers For each of these four wind conditions carry out the following computations Find the 17 trajectories whose initial angles are multiples of 5 degrees that is 90 Jew36 radians k 1 21 1 Plot all 17 trajectories on one gure Determine which of these trajectories has the greatest downrange distance For that trajectory report the initial angle in degrees the ight time the downrange distance the impact velocity and the number of steps required by the ODE solveri Which of the four wind conditions requires the most computation Why 1 1n the 1968 Olympic games in Mexico City Bob Beamon established a world record with a longjump of 8190 meters This was 0180 meters longer than the previous world record Since 1968 Beamon s jump has been exceeded only once in competition by Mike Powell s jump of 895 meters in Tokyo in 1991 After Beamon s remarkable jump some people suggested that the lower air resistance at Mexico Cityls 7400 ft altitude was a contributing factor This problem examines that possibility Exercises 41 The mathematical model is the same as the cannonball trajectory in the previous exercise The xed Cartesian coordinate system has a horizontal zaxis a vertical y axis and an origin at the takeoff board The jumper7s initial velocity has magnitude v0 and makes an angle With respect to the zaxis of 90 radiansi The only forces acting after takeoff are gravity and the aerodynamic drag D Which is proportional to the square of the magnitude of the velocity There is no wind The equations describing the jumpers motion are ivcost9 y39vsint9 D 67gCOSO 137iigsint9 v m The drag is Constants for this exercise are the acceleration of gravity 9 981 ms2 the mass m 80 kg the drag coefficient 5 072 the jumpers crosssectional area 8 050 m2 and the takeoff angle 90 2250 7r8 radiansi Compute four different jumps With different values for initial velocity v0 and air density p The length of each jump is 1tf Where the air time tf is determined by the condition ytf 0 a Nominal jump at high altitude v0 10 ms and p 0 94 kgmgi b Nominal jump at sea level v0 10 ms and p 1 29 kgmgi c Sprinter7s approach at high altitude p 0 94 kgmgi Determine v0 so that the length of the jump is Beamon s record 890 m d Sprinter s approach at sea level p 1 29 kgmg and v0 is the value determined in Present your results by completing the following table v0 thetaO rho distance 100000 225000 09400 100000 225000 12900 225000 09400 8 9000 22 5000 12900 Which is the more important the air density or the jumpers initial velocity A pendulum is a stiff bar of length L supported at one end by a frictionless pini lf gravity is the only force acting on the pendulum its oscillation is described by the ODE 5 7gL sint9 Here 9 is the angular position of the bar With 9 0 if the bar is hanging down from the pin and t9 7r if the bar is precariously balanced above the pin Take L 12 inches and g 39609 inchess2i If the pendulum is released 42 Chapter 7 Ordinary Differential Equations from rest with an initial angle 90 the period T of the oscillatory motion is given by Two 4L912K902 where is the complete elliptic integral of the rst kind given by 7r2 u d v 0 l 7 sin2 u sin2 1 a Use dsolve in the Symbolic Toolbox to solve the differential equation for small oscillations small 9 and therefore show that Two approaches 27rLg12 as 90 approaches zero For parts b through e do not make the assumption that 9 is small and so do not replace sini9 by 9 b Explain on physical grounds why Two approaches in nity as 90 ap proaches 7r c Con rm numerically that the analytical formula given above for the period is correct by writing a function that uses QUAD4 and computes the period Two for 0 S 90 S 0 99997ri Make a plot of this function Warning you have to be careful if 90 is near 7r d Compute the motion over one period for several different values of 90 including values near 0 and near 7r Superimpose the phase plane plots of the solutions for several different values of 90 on one graphi What is the effect of burning fossil fuels on the carbon dioxide in the earthls atmosphere Even though today C02 accounts for only about 350 parts per million of the atmosphere any increase has profound implications for our climate An informative background article is available at a Web site maintained by the Lighthouse Foundation A model developed by Jr C Gr Walker 9 was brought to our attention by Eric Rodeni The model simulates the interaction of the various forms of carbon that are stored in three regimes the atmosphere the shallow ocean and the deep ocean The ve principal variables in the model are all functions of time p partial pressure of carbon dioxide in the atmosphere as total dissolved carbon concentration in the shallow ocean 00 total dissolved carbon concentration in the deep ocean as alkalinity in the shallow ocean ad alkalinity in the deep ocean Three additional quantities are involved in equilibrium equations in the shal low oceani hs hydrogen carbonate in the shallow ocean cs carbonate in the shallow ocean p5 partial pressure of gaseous carbon dioxide in the shallow oceani Exercises 43 The rate of change of the ve principal variables is given by ve ordinary differential equations The exchange between the atmosphere and the shallow ocean involves a constant characteristic transfer time d and a source term f t 7 P5 P M dt 7 d 1 The equations describing the exchange between the shallow and deep oceans involve v5 and vd the volumes of the two regimesi 7w 7 7 h 7 w vid k1 7 0390 7 MW dis vi ad 7 a5w 7 k2 vid kg 7 ad 7 10w The equilibrium between carbon dioxide and the carbonates dissolved in the shallow ocean is described by three nonlinear algebraic equations h 7 as 7 a 7 kgaseas 7 s 7 M C 7 as 7 hs 5 7 2 h2 P5 44 55 The numerical values of the constants involved in the model are d 7 864 m 7 495102 2 7 495102 US 7 012 W 7 123 w 7 103 k1 7 21910 4 kg 7 61210 5 1 7 0997148 M 7 67910 2 44 Chapter 7 Ordinary Differential Equations The source term describes the burning of fossil fuels in the modern industrial erai We will use a time interval that starts about a thousand years ago and extends a few thousand years into the future 1000 S t S 5000 The initial values at t 1000 p 1100 a 2101 ad 2123 a 2120 ad 2 26 represent preindustrial equilibrium and remain nearly constant as long as the source term is zero The following table describes one scenario for a source term that models the release of carbon dioxide from burning fossil fuels especially gasoline The amounts begin to be signi cant after 1850 peak near the end of this century and then decrease until the supply is exhausted year rate 1000 0 1850 0 0 1950 1 0 1980 4 0 2000 5 0 2050 8 0 2080 10 0 2100 10 5 2120 10 0 2150 8 0 2225 3 5 2300 2 0 2500 0 0 5000 00 Figure 7110 shows this source term and its effect on the atmosphere and the ocean The three graphs in the lower half of the gure show the atmospheric shallow ocean and deep ocean carbon The two alkalinity values are not plotted at all because they are almost constant throughout this entire simu lation Initially the carbon in the three regimes is nearly at equilibrium and so the amounts hardly change before 1850 Over the period 1850 S t S 2500 the upper half of gure 7110 shows the additional carbon produced by burning fossil fuels entering the system and the lower half shows the system response The atmosphere is the rst to be affected showing more than a fourfold increase in 500 years Almost half of the carbon is then slowly transferred to the shallow ocean and eventually to the deep ocean Exercises 45 Carbon in the atmosphere and ocean 15 l l l l l 10 i E a 5 fossil fuel l l l l l 1000 1500 2000 2500 3000 3500 4000 4500 5000 7 H carbon shallow 1500 2000 2500 3000 3500 4000 4500 5000 time yr Figure 710 Carbon in the atmosphere and ocean a Reproduce gure 710 Use pchiptx to interpolate the fuel table and ode23tx With the default tolerances to solve the differential equations b How do the amounts of carbon in the three regimes at year 5000 compare With the amounts at year 1000 c When does the atmospheric C02 reach its maximum d These equations are mildly stiff because the various chemical reactions take place on very different time scales If you zoom in some portions of the graphs you should see a characteristic sawtooth behavior caused by the small time steps required by ode23txl Find such a region e Experiment With other MATLAB ODE solvers including ode23 ode45 ode113 ode235 and odelEsl Try various tolerances and report computa tional costs by using something like odeset Re1Tol 1e6 AbsTol 1e6 stats on Which method is preferable for this problem i This problem makes use of quadrature ordinary differential equations and zero nding to study a nonlinear boundaryvalue probleml The function 46 Chapter 7 Ordinary Differential Equations is de ned on the interval 0 S I S l by y y2 71 90 0 91 1 This problem can be solved four different ways Plot the four solutions ob tained on a single figure7 using subplot 2 2 1 subplot 2 2 4 a Shooting methodi Suppose we know the value of 77 yO Then we could use an ODE solver like ode23tx or ode45 to solve the initial value problem yn y2 71 90 0 y 0 77 on the interval 0 S I S 1 Each value of 77 determines a different solution yz 77 and corresponding value for yl 77 The desired boundary condition yl 1 leads to the definition of a function of 77 f0 y177 1 Write a MATLAB function whose argument is 77 This function should solve the ODE initial problem and return Then use fzero or fzerotx to find a value 77 so that 0i Finally7 use this 77 in the initial value problem to get the desired Report the value of 77 you obtain b Quadrature Observe that y y2 7 1 can be written iltlty gti ygto dz 2 3 This means that the expression 7 9 y3 H 7 2 7 3 y is actually constanti Because y0 07 we have y 0 V217 So if we could find the constant H the boundaryvalue problem would be converted into an initial value problemi Integrating the equation dz 1 2H m gives y I hm dy 0 where 1 hm 3 2Hy 3y Exercises 47 This together with the boundary condition yl 1 leads to the de nition of a function 9H 9HA Mali dy 7 1 You need two MATLAB functions one that computes hyH and one that computes gHi They can be two separate M les but a better idea is to make hyH an inline function within 9Hi The function 9H should use quadtx to evaluate the integral of hy H The parameter H is passed as an extra argument from 9 through quadtx to h Then fzerotx can be used to nd a value li so that 90 0 Finally this li provides the second initial value necessary for an ODE solver to compute Report the value of li you obtain c and d Nonlinear nite differences Partition the interval into n 1 equal subintervals with spacing h ln 1 z ih i0iunl Replace the differential equation by a nonlinear system of difference equations involving n unknowns yl yg i i i yn yi1 Zyi l yiil h2yi2717 i17n The boundary conditions are yo 0 and yn1 1 A convenient way to compute the vector of second differences involves the nbyn tridiagonal matrix A with 727s on the diagonal 17s on the super and subdiagonal and 07s elsewhere You can generate a sparse form of this matrix with D onesn1 spdiagse 2e e 1 0 1nn p The boundary conditions yo 0 and yn1 1 can be represented by the nvector b with bi 0i l i i i n 7 l and In 1 The vector formulation of the nonlinear difference equation is Aybh2lty271gt where y2 is the vector containing the squares of the elements of y that is the MATLAB elementbyelement power y quot2 There are at least two ways to solve this systemi c Linear iterationi This is based on writing the difference equation in the form Ay h2y2 70712 Start with an initial guess for the solution vector y The iteration consists of plugging the current y into the righthand side of this equation and then solving the resulting linear system for a new y This makes repeated use of the sparse backslash operator with the iterated assignment statement 48 Chapter 7 Ordinary Differential Equations y Ah 2y 2 1 b It turns out this iteration converges linearly and provides a robust method for solving the nonlinear difference equations Report the value of n you use and the number of iterations require d Newtonls method This is based on writing the difference equation in the form Fy Ayb7h2y21 0 Newtonls method for solving Fy 0 requires a manyvariable analog of the derivative F The analog is the Jacobian the matrix of partial derivatives 7 BF Byj J A 7 h2diag2y ln MATLAB one step of Newtonls method would be F My b h 2y 2 1 A h 2spdiags2y0nn y y JF lt With a good starting guess Newtonls method converges in a handful of iter ationsi Report the value of n you use and the number of iterations requiredi The double pendulum is a classical physics model system that exhibits chaotic motion if the initial angles are large enough The model shown in gure 7 11 involves two weights or 70175 attached by weightless rigid rods to each other and to a xed pivoti There is no friction so once initiated the motion continues forever The motion is fully described by the two angles 91 and 92 that the rods make with the negative y axis Figure 711 Double pendulum Exercises 49 Let m1 and 7212 be the masses of the bobs and Z1 and 2 be the lengths of the rods The positions of the bobs are 11 Z1 sin91 yl 7Z1 cos 91 I2 Z1 sin91 Z2 sin92 y2 7Z1 cos91 7 2 cos 92 The only external force is gravity denoted by g Analysis based on the La grangian formulation of classical mechanics leads to a pair of coupled second order nonlinear ordinary differential equations for the two angles 91t and 92 m1 m2 191 m2 292 cos 91 7 62 9m1 m2 Sir191 7712529 Sin 91 i 92 m2Z191 cos 91 7 92 m2Z292 7gm2 sin 92 721259 sin 91 7 92 To rewrite these equations as a rst order system introduce the 4byl column vector ut u 1917 1927 1917 92lT With m1 m2 Z1 Z2 1 c cos ul 7 u2 and s sin ul 7 u2 the equations become ul usy 13912 U47 2ugcu4 7gsinu17sui cugu4 7gsinu2su Let M denote the 4by4 mass matrix 1000 0100 M002c 0051 and let f denote the 4 byl nonlinear force function us u4 7y sin ul 7 sui 7y sin u2 311 In matrixvector notation the equations are simply Muf 50 Chapter 7 Ordinary Differential Equations This is an implicit system of differential equations involving a nonconstant nonlinear mass matrix The double pendulum problem is usually formulated without the mass matrix but larger problems with more degrees of freedom are frequently in implicit formi In some situations the mass matrix is singular and it is not possible to write the equations in explicit formi The NCM M le swinger provides an interactive graphical implemention of these equations The initial position is determined by specifying the starting coordinates of the second bob 12 yg either as arguments to swinger or by using the mouse In most situations this does not uniquely determine the starting position of the rst bob but there are only two possibilities and one of them is chosen arbitrarily The initial velocities 91 and 392 are zero The numerical solution is carried out by ode23 because our textbook code ode23tx cannot not handle implicit equations The call to ode23 involves using odeset to specify the functions that generate the mass matrix and do the plotting opts odeset mass swingmass outputfcn swingplot ode23 swingrhs tspan uO opts The mass matrix function is function M swingmasstu c cosu1u2 M10000100002c00c1 The driving force function is function f swingrhstu g 1 sinu1u2 Ham It would be possible to have just one ode function that returns Mf but we want to emphasize the implicit facilityi An internal function swinginit converts a speci ed starting point zy to a pair of angles 91 92 If I y is outside the circle i gthb then the pendulum cannot reach the speci ed point In this case we straighten out the pendulum with 91 92 and point it in the given direction If I y is inside the circle of radius 2 we return one of the two possible con gurations that reach to that point Here are some questions to guide your investigation of swingeri a When the initial point is outside the circle of radius 2 the two rods start out as one If the initial angle is not too large the double pendulum continues to act pretty much like a single pendulumi But if the initial angles are large enough chaotic motion insuesi Roughly what initial angles lead to chaotic motioni b The default initial condition is u3 u4 2gsinu1su4 2 gsinu2su3 2 Exercises 51 swinger O 862 0 994 Why is this orbit interesting Can you nd any similar orbits c Run swinger for a While then click on its stop button Go to the MATLAB command line and type get gcf userdata 1 What is returned d Modify swinginit so that When the initial point is inside the circle of radius 2 the other possible initial con guration is chosen e Modify swinger so that masses other than m1 m2 1 are possible f Modify swinger so that lengths other than 1 Z2 1 are possible This is tricker than changing the masses because the initial geometry is involved g What role does gravity play How would the behavior of a double pendu lum change if you could take it to the Moon How does changing the value of g in swingrhs affect the speed of the graphics display the step sizes chosen by the ode solver and the computed values of t h Combine swingmass and swingrhs into one function swingodei Elimi nate the mass option and use ode23tx instead of ode231 i Are these equations stiff This is a difficult question The statement swinger02 tries to deli cately balance the pendulum above its pivot point The pendulum does stay there for a While but then loses its balance Observe the value oft displayed in the title for swinger 02 What force knocks the pendulum away from the vertical position At What value oft does this force become noticeable 52 Chapter 7 Ordinary Differential Equations 11 E Bibliography Ui Mi ASCHER AND L RY PETZOLD Computer Methods for Ordinary Di er ential Equations and Di erentialAlqebraic Equations SIAM 1998 KEY BRENAN SlLi CAMPBELL AND LiRl PETZOLD The Numerical Solution of Initial Value Problems in Di erentialAlqebraic Equations SIAM 1996i Pi BOGACKI AND L F SHAMPINE A pair of RunqeKutta formulas Appl Math Letters 2 1989 pp 179 R Ml CORLESS G H GONNET D E G HARE D J JEFFREY AND D E KNUTH On the Lambert W Function Advances in Computational Mathematics volume 5 1996 pp 3297359 http wwwapmaths uwo ca rcorlessframesPAPERSLambertw LIGHTHOUSE FOUNDATION http wwwlighthousefoundationorglighthousefoundationorg engexplorerartike100294enghtml L F SHAMPINE Numerical solution of ordinary di erential equations Chap man ancl Hall New York 1994 Lo Fi SHAMPINE AND Ml WY REICHELT The MATLAB ODE Suite SIAM J Scienti c Computing 18 1997 pp 1722 C SPARROW The Lorenz Equations Bifurcations Chaos and Strange Attrac tors SpringerVerlag New York 1982 269 pp J C Go WALKER Numerical adventures with geochemical cycles Oxford Uni versity Press New York 1991


Buy Material

Are you sure you want to buy this material for

25 Karma

Buy Material

BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.


You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

Why people love StudySoup

Steve Martinelli UC Los Angeles

"There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

Janice Dongeun University of Washington

"I used the money I made selling my notes & study guides to pay for spring break in Olympia, Washington...which was Sweet!"

Bentley McCaw University of Florida

"I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

Parker Thompson 500 Startups

"It's a great way for students to improve their educational experience and it seemed like a product that everybody wants, so all the people participating are winning."

Become an Elite Notetaker and start selling your notes online!

Refund Policy


All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email


StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here:

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.