[Solved] MATLAB sessions: Laboratory 4 MAT 275 Laboratory 4

MATLAB sessions: Laboratory 4MAT 275 Laboratory 4MATLAB solvers for First-Order IVPIn this laboratory session we will learn how to1. Use MATLAB solvers for solving scalar IVP2. Use MATLAB solvers for solving higher order ODEs and systems of ODES.First-Order Scalar IVPConsider the IVP

y0 = t y,y(0) = 1.(L4.1)The exact solution is y(t) = t 1 + 2et. A numerical solution can be obtained using various MATLABsolvers. The standard MATLAB ODE solver is ode45. The function ode45 implements 4/5th orderRunge-Kutta method. Type help ode45 to learn more about it.Basic ode45 UsageThe basic usage of ode45 requires a function (the right-hand side of the ODE), a time interval on whichto solve the IVP, and an initial condition. For scalar first-order ODEs the function may often be specifiedusing the inline MATLAB command. A complete MATLAB solution would read:1 f = inline(t-y,t,y);2 [t,y] = ode45(f,[0,3],1);3 plot(t,y)(line numbers are not part of the commands!) Line 1 defines the function f as a function of t and y, i.e., f(t, y) = t y. This is the right-handside of the ODE (L4.1). Line 2 solves the IVP numerically using the ode45 solver. The first argument is the function f,the second one determines the time interval on which to solve the IVP in the form[initial time, final time], and the last one specifies the initial value of y. The output of ode45consists of two arrays: an array t of discrete times at which the solution has been approximated,and an array y with the corresponding values of y. These values can be listed in the CommandWindow as[t,y]ans =0 1.00000.0502 0.95220.1005 0.90930.1507 0.87090.2010 0.8369..2.9010 2.01092.9257 2.03302.9505 2.05512.9752 2.07733.0000 2.0996c 2016 Stefania Tracogna, SoMSS, ASU 1MATLAB sessions: Laboratory 4Since the output is quite long we printed only some selected values.For example the approximate solution at t 2.9257 is y 2.0330. Unless specific values of y areneeded it is better in practice to simply plot the solution to get a sense of the behavior of thesolution. Line 3 thus plots y as a function of t in a figure window. The plot is shown in Figure L4a.0 0.5 1 1.5 2 2.5 30.60.811.21.41.61.822.2Figure L4a: Solution of (L4.1).Error Plot, Improving the AccuracyError plots are commonly used to show the accuracy in the numerical solution. Here the error is thedifference between the exact solution y(t) = t1+ 2et and the numerical approximation obtained fromode45. Since this approximation is only given at specified time values (contained in the array t) we onlyevaluate this error at these values of t:err = t-1+2*exp(-t)-yerr =1.0e-005 *00.02780.04070.0162-0.0042-0.0329-0.0321-0.0313-0.0305-0.0298(in practice the exact solution is unknown and this error is estimated, for example by comparing thesolutions obtained by different methods). Again, since the error vector is quite long we printed only afew selected values. Note the 1.0e-005 at the top of the error column. This means that each componentof the vector err is less than 105in absolute value.A plot of err versus t is more revealing. To do this note that errors are usually small so it is bestto use a logarithmic scale in the direction corresponding to err in the plot. To avoid problems withnegative numbers we plot the absolute value of the error (values equal to 0, e.g. at the initial time, arenot plotted):semilogy(t,abs(err)); grid on;c 2016 Stefania Tracogna, SoMSS, ASU 2MATLAB sessions: Laboratory 40 0.5 1 1.5 2 2.5 3108107106105Figure L4b: Error in the solution of (L4.1) computed by ode45.See Figure L4b. Note that the error level is about 106. It is sometimes important to reset the defaultaccuracy ode45 uses to determine the approximation. To do this use the MATLAB odeset commandprior to calling ode45, and include the result in the list of arguments of ode45:f = inline(t-y,t,y);options = odeset(RelTol,1e-10,AbsTol,1e-10);[t,y] = ode45(f,[0,3],1,options);err = t-1+2*exp(-t)-y;semilogy(t,abs(err))See Figure L4c.0 0.5 1 1.5 2 2.5 31016101510141013101210111010Figure L4c: Error in the solution of (L4.1) computed by ode45 with a better accuracy.Parameter-Dependent ODEWhen the function defining the ODE is complicated, rather than entering it as inline, it is moreconvenient to enter it as a separate function file, included in the same file as the calling sequence ofode45 (as below) or saved in a separate m-file.In this Section we will look at an example where the ODE depends on a parameter.Consider the IVP

y0 = a (y et) et,y(0) = 1.(L4.2)with exact solution y(t) = et(independent of the parameter a!). An implementation of the MATLABsolution in the interval [0, 3] follows.c 2016 Stefania Tracogna, SoMSS, ASU 3MATLAB sessions: Laboratory 41 function ex_with_param2 t0 = 0; tf = 3; y0 = 1;3 a = 1;4 [t,y] = ode45(@f,[t0,tf],y0,[],a);5 disp([y( num2str(t(end)) ) = num2str(y(end))])6 disp([length of y = num2str(length(y))])7 %-8 function dydt = f(t,y,a)9 dydt = -a*(y-exp(-t))-exp(-t); Line 1 must start with function, since the file contains at least two functions (a driver + afunction). Line 2 sets the initial data and the final time. Line 3 sets a particular value for the parameter a. In line 4 the parameter is passed to ode45 as the 5th argument (the 4th argument is reserved forsetting options such as the accuracy using odeset, see page 3, and the placeholder [] must beused if default options are used).Correspondingly the function f defined in lines 8-9 must include a 3rd argument corresponding tothe value of the parameter. Note the @f in the argument of ode45. This is necessary when thefunction is not defined as inline. See the help on ode45 for more information. On line 5 the value of y(3) computed by ode45 is then displayed in a somewhat fancier form thanthe one obtained by simply entering y(end). The command num2string converts a number to astring so that it can be displayed by the disp command.The m-file ex with param.m is executed by entering ex with param at the MATLAB prompt. Theoutput is>> ex_with_paramy(3) = 0.049787length of y = 45 The additional line 6 in the file lists the length of the array y computed by ode45. It is interestingto check the size of y obtained for larger values of a. For example for a = 1000 we obtain>> ex_with_paramy(3) = 0.049792length of y = 3621This means that ode45 needed to take smaller step sizes to cover the same time interval comparedto the case a = 1, even though the exact solution is the same!Not all problems with a common solution are the same! Some are easier to solve than others.When a is large the ODE in (L4.2) is said to be stiff. Stiffness has to do with how fast nearbysolutions approach the solution of (L4.2), see Figure L4d.F Other MATLAB ODE solvers are designed to better handle stiff problems. For example replace ode45with ode15s in ex with param.m (without changing anything else) and set a = 1000:4 [t,y] = ode15s(@f,[t0,tf],y0,[],a);>> ex_with_paramy(3) = 0.049787length of y = 18c 2016 Stefania Tracogna, SoMSS, ASU 4MATLAB sessions: Laboratory 4Figure L4d: Direction field and sample solutions in the t-y window [0, 0.1] [0.9, 1] as obtained usingDFIELD8: a = 1 (left) and a = 1000 (right).A solution exhibiting blow up in finite timeConsider the Differential Equationdydt = 1 + y2, y(0) = 0The exact solution of this IVP is y = tan t and the domain of validity is [0,2). Lets see what happenswhen we try to implement this IVP using ode45 in the interval [0, 3].>> f=inline(1+y^2,t,y);>> [t,y]=ode45(f,[0,3],0);Warning: Failure at t=1.570781e+000. Unable to meet integrationtolerances without reducing the step size below the smallest valueallowed (3.552714e-015) at time t.> In ode45 at 371The MATLAB ode solver gives a warning message when the value of t = 1.570781 is reached. This isextremely close to the value of /2 where the vertical asymptote is located.If we enter plot(t,y) we obtain Figure L4e on the left (note the scale on the y-axis), however, if we usexlim([0,1.5]), we can recognize the graph of y = tan t.0 0.5 1 1.5 2024681012141618x 10130 0.5 1 1.5051015Figure L4e: Solution of y0 = 1 + y2, y(0) = 0 without restrictions on the axis, and with xlim([0,1.5])c 2016 Stefania Tracogna, SoMSS, ASU 5MATLAB sessions: Laboratory 4Higher-Order and Systems of IVPsWe show here how to extend the use of ode45 to systems of first-order ODEs (the same holds for othersolvers such as ode15s). Higher-order ODEs can first be transformed into a system of first-order ODEsto fit into this framework. We will see later how to do this.As an example consider the predator-prey system (Lotka-Volterra) representing the evolution of twopopulations. u1 = u1(t) and u2 = u2(t):du1dt = au1 bu1u2,du2dt = cu2 + du1u2(L4.3)with initial populations u1(0) = 10 and u2(0) = 60. The parameters a, b, c, and d are set to a = 0.8,b = 0.01, c = 0.6, and d = 0.1. The time unit depends on the type of populations considered.Although the ODE problem is now defined with two equations, the MATLAB implementation isvery similar to the case of a single ODE, except that vectors must now be used to describe the unknownfunctions.1 function ex_with_2eqs2 t0 = 0; tf = 20; y0 = [10;60];3 a = .8; b = .01; c = .6; d = .1;4 [t,y] = ode45(@f,[t0,tf],y0,[],a,b,c,d);5 u1 = y(:,1); u2 = y(:,2); % y in output has 2 columns corresponding to u1 and u26 figure(1);7 subplot(2,1,1); plot(t,u1,b-+); ylabel(u1);8 subplot(2,1,2); plot(t,u2,ro-); ylabel(u2);9 figure(2)10 plot(u1,u2); axis square; xlabel(u_1); ylabel(u_2); % plot the phase plot11 %-12 function dydt = f(t,y,a,b,c,d)13 u1 = y(1); u2 = y(2);14 dydt = [ a*u1-b*u1*u2 ; -c*u2+d*u1*u2 ]; In line 2 the 2 1 vector y0 defines the initial condition for both u1 and u2. In line 4 the parameters are passed to the ODE solver ode45 as extra arguments (starting fromthe 5th), as many as there are parameters in the problem (4 here).The output array y of ode45 now has 2 columns, corresponding to approximations for u1 and u2,respectively, instead of a single one. In line 5 these quantities are therefore retrieved and stored in arrays u1 and u2, which are descriptivenames. The part of the program defining the ODE system includes lines 11-14. Note that all the parametersappearing as arguments of ode45 must appear as arguments of the function f. For a specific valueof t the input y to f is a 2 1 vector, whose coefficients are the values of u1 and u2 at time t.Rather than referring to y(1) and y(2) in the definition of the equations on line 14, it is bestagain to use variable names which are easier to identify, e.g., u1 and u2. Line 14 defines the right-hand sides of the ODE system as a 2 1 vector: the first coefficient isthe first right-hand side ( du1dt ) and the second coefficient the second right-hand side ( du2dt ). Lines 6-10 correspond to the visualization of the results. To plot the time series of u1 and u2, wecreate a 2 1 array of subplots. Because the scales of u1 and u2 are different, it is best using twodifferent graphs for u1 and u2 here. Type help subplot to learn more about it. On a differentfigure, we then plot the phase plot representing the evolution of u2 in terms of u1. Note that u1and u2 vary cyclically. The periodic evolution of the two populations becomes clear from the closedcurve u2 vs. u1 in the phase plot.c 2016 Stefania Tracogna, SoMSS, ASU 6MATLAB sessions: Laboratory 40 5 10 15 2024681012u10 5 10 15 20406080100120140u22 4 6 8 10 12405060708090100110120130140u1u2Figure L4f: Lotka-Volterra example.Reducing a Higher-Order ODENumerical solution to IVPs involving higher order ODEs homogeneous or not, linear or not, can beobtained using the same MATLAB commands as in the first-order by rewriting the ODE in the form ofa system of first order ODEs.Lets start with an example. Consider the IVPd2ydt2+ 7dydt + 5y = sin t, with y(0) = 0.5,dydt (0) = 0.5. (L4.4)To reduce the order of the ODE we introduce the intermediate unknown function v =dydt . As a resultdvdt =d2ydt2 so that the ODE can be written dvdt + 7v + 5y = sin t. This equation only involves first-orderderivatives, but we now have two unknown functions y = y(t) and v = v(t) with two ODEs. ForMATLAB implementations it is necessary to write these ODEs in the form ddt = . . .. Thusd2ydt2+ 7dydt + 5y = sin t ( dydt = v,dvdt = sin t 7v 5y.(L4.5)Initial conditions from (L4.4) must also be transformed into initial conditions for y and v. Simply,y(0) = 0.5,dydt (0) = 0.5 (y(0) = 0.5,v(0) = 0.5.(L4.6)EXERCISESInstructions: For your lab write-up follow the instructions of LAB 1.1. (a) Modify the function ex with 2eqs to solve the IVP (L4.4) for 0 t 50 using the MATLABroutine ode45. Call the new function LAB04ex1.Let [t,Y] (note the upper case Y) be the output of ode45 and y and v the unknown functions.Use the following commands to define the ODE:function dYdt= f(t,Y)y=Y(1); v=Y(2);dYdt = [v; sin(t)-7*v-5*y];Plot y(t) and v(t) in the same window (do not use subplot), and the phase plot showing vvs y in a separate window.Add a legend to the first plot. (Note: to display v(t) = y0(t), use v(t)=y(t)).Add a grid. Use the command ylim([-1.2,1.2]) to adjust the y-limits for both plots.Adjust the x-limits in the phase plot so as to reproduce the pictures in Figure L4g.Include the M-file in your report.c 2016 Stefania Tracogna, SoMSS, ASU 7MATLAB sessions: Laboratory 40 5 10 15 20 25 30 35 40 45 50-1-0.500.51y(t)v(t)-0.5 0 0.5 1y-1-0.500.51v=yFigure L4g: Time series y = y(t) and v = v(t) = y0(t) (left), and phase plot v = y0 vs. y for (L4.4).(b) For what (approximate) value(s) of t does y reach a local maximum in the window 0 t 50?Check by reading the matrix Y and the vector t. Note that, because the M-file LAB04ex1.mis a function file, all the variables are local and thus not available on the Command Window.To read the matrix Y and the vector t, you need to modify the M-file by adding the line[t Y].Do not include the whole output in your lab write-up. Include only the values necessary toanswer the question.(c) What seems to be the long term behavior of y?(d) Modify the initial conditions to y(0) = 2, v(0) = 3 and run the file LAB04ex1.m with themodified initial conditions. Based on the new graphs, determine whether the long termbehavior of the solution changes. Explain. Include the pictures with the modified initialconditions to support your answer.Nonlinear ProblemsNonlinear problems do not present any additional difficulty from an implementation point of view (theymay present new numerical challenges for integration routines like ode45).EXERCISES2. (a) Consider the modified problemd2ydt2+ 7y2dydt + 5y = sin t, with y(0) = 0.5,dydt (0) = 0.5. (L4.7)The ODE (L4.7) is very similar to (L4.4) except for the y2term in the left-hand side. Becauseof the factor y2the ODE (L4.7) is nonlinear, while (L4.4) is linear. There is however verylittle to change in the implementation of (L4.4) to solve (L4.7). In fact, the only thing thatneeds to be modified is the ODE definition.Modify the function defining the ODE in LAB04ex1.m. Call the revised file LAB04ex2.m. Thenew function M-file should reproduce the pictures in Fig L4h.Include in your report the changes you made to LAB04ex1.m to obtain LAB04ex2.m.(b) Compare the output of Figs L4g and L4h. Describe the changes in the behavior of the solutionin the short term.(c) Compare the long term behavior of both problems (L4.4) and (L4.7), in particular the amplitudeof oscillations.(d) Modify LAB04ex2.m so that it solves (L4.7) using Eulers method with N = 1000 in theinterval 0 t 50 (use the file euler.m from LAB 3 to implement Eulers method; donot delete the lines that implement ode45). Let [te,Ye] be the output of euler, and notec 2016 Stefania Tracogna, SoMSS, ASU 8MATLAB sessions: Laboratory 40 5 10 15 20 25 30 35 40 45 50-1-0.500.51y(t)v(t)-0.5 0 0.5 1y-1-0.500.51v=yFigure L4h: Time series y = y(t) and v = v(t) = y0(t) (left), and phase plot v = y0 vs. y for (L4.7).that Ye is a matrix with two columns from which the Eulers approximation to y(t) must beextracted. Plot the approximation to the solution y(t) computed by ode45 (in black) and theapproximation computed by euler (in red) in the same window (you do not need to plot v(t)nor the phase plot). Are the solutions identical? Comment. What happens if we increase thevalue of N?Include the modified M-file in your report.3. Solve numerically the IVPd2ydt2+ 7ydydt + 5y = sin t, with y(0) = 0.5,dydt (0) = 0.5in the interval 0 t 50. Include the M-file in your report.Is the behavior of the solution significantly different from that of the solution of (L4.7)?Is MATLAB giving any warning message? Comment.A Third-Order ProblemConsider the third-order IVPd3ydt3+ 7y2d2ydt2+ 14y

dydt 2+ 5dydt = cost, with y(0) = 0.5,dydt (0) = 0.5,d2ydt2(0) = 1.625. (L4.8)Introducing v =dydt and w =dy2dt2 we obtain dvdt = w and dwdt =d3ydt3 = cost7y2w14yv2 5v. Moreover,v(0) = dydt (0) = 0.5 and w(0) = d2ydt2 (0) = 1.625. Thus (L4.8) is equivalent todydt = v,dvdt = w,dwdt = cost 7y2w 14yv2 5vwithy(0) = 0.5,v(0) = 0.5,w(0) = 1.625.(L4.9)4. (a) Write a function M-file that implements (L4.8) in the interval 0 t 50. Note that theinitial condition must now be in the form [y0,v0,w0] and the matrix Y, output of ode45,has now three columns (from which y, v and w must be extracted). On the same figure, plotthe three time series and, on a separate window, plot the phase plot usingfigure(2); plot3(y,v,w,k.-);grid on; view([-40,60])xlabel(y); ylabel(v=y); zlabel(w=y);c 2016 Stefania Tracogna, SoMSS, ASU 9MATLAB sessions: Laboratory 40 5 10 15 20 25 30 35 40 45 50-1-0.500.51y(t)v(t)w(t)-41.50.6 10.40.5 0.2v=y y0 0-2-0.2-0.5-0.4-1 -0.6w=y02Figure L4i: Time series y = y(t), v = v(t) = y0(t), and w = w(t) = y0(t) (left), and 3D phase plot v = y0vs. y vs w = y00 for (L4.8) (rotated with view = [-40,60]).Do not forget to modify the function defining the ODE.The output is shown in Fig. L4i. The limits in the vertical axis of the plot on the left weredeliberately set to the same ones as in Fig. L4h for comparison purposes using the MATLABcommand ylim([-1.2,1.2]).Include the M-file in you report.(b) Compare the output of Figs L4h and L4i. What do you notice?(c) Differentiate the ODE in (L4.7) and compare to the ODE in (L4.8).(d) Explain why the solution of (L4.7) also satisfies the initial conditions in (L4.8). Hint: Substitutet = 0 into (L4.7) and use the initial conditions for y and v.c 2016 Stefania Tracogna, SoMSS, ASU 10