Numerical Applications to Differential Equations
Numerical Applications to Differential Equations MA 302
Popular in Course
Popular in Mathematics (M)
verified elite notetaker
This 99 page Class Notes was uploaded by Braeden Lind on Thursday October 15, 2015. The Class Notes belongs to MA 302 at North Carolina State University taught by Robert White in Fall. Since its upload, it has received 12 views. For similar materials see /class/223678/ma-302-north-carolina-state-university in Mathematics (M) at North Carolina State University.
Reviews for Numerical Applications to Differential Equations
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
Differential Equation Model Continuous Version of Newton s Law of Cooling y Cysur y and y0 yo Discrete Version of Newton s Law of Cooling yil yi cysur yidt We can solve for yil yi1 yi 0ysur gt39iquot dt 1 cdtyi cysurdt Method of Solution Find c by using two observations y0 200 and y1 195 Now apply the model with i O to get y01 1 cdty0 cysurdt 195 1 c 1200 c701 So solve for c to get c 195 20070 200 5130 126 Experiment with different 11 KK 481632 Compare the errors in the exact and numerical solutions Exact solution at t 50 is y50 8900035 Numerical solution at t 50 is Ynl Error 890035 Ynl 955137 458938 224524 lllO24 forn4 forn8 forn 16 for n 32 Matlab Implementation Use function le feulm and m le eulercm function feul feultx feul 1070 x your name your student number lesson number clear y1 200 T 50 KK 100 h TKK t1 0 for k lzKK yk1 yk hfeultk Yk tk1 tk h end plottY title39your name your student number lesson number39 xlabel39time39 ylabel39temperature39 temperature 160 your name your student number lesson number Differential Equation Model Let S 7727 W 161 q 837510A6 y10 4y20 11 and y30 4 y139 sY2 39 y1y2 y1 q 3 12 y239 y2 y1 y2 y3s and y339 Wy1 39 Y3 Matlab39s 0de23s and 0de15s can be used to solve such systems The steady state solutions which are defined by all three derivatives being set equal to zero so that 0 sY2 39 y1y2 y1 39 q 3 12 f1y1a Y2 Y3 a 0 39 Y2 39 y1 Y2 Y35 f2y1a Y2 Y3 and 0 WYI 39 Y3 f3 h Y2 Y3 The three positive steady state solutions are y1 y3 48818 and y 99796 These may be found by using the symbolic Matlab command solve as follows EDU y1y2solve39y2 y1y2y1 qy1y103939y2y1y2y1039 The Jacobian of f1 f2 f339 is n f1 f1 4 1612y1 50 h 0 1ng f2 f2 y2s 1 yls ls J2 Jab y w 0 w The eigenvalues of J evaluated at the steady state solutions are easily found by using the Matlab command eig The following code is in the Matlab m le called 0regjacm y1 48818 y2 99796 Y3 yl s 7727 w 161 q 83751OA 6 J s y21 q2y1 s1 y1 O y2s lyls 15 W 0 w eigJ EDU oregjac ans 257146 187473 00013 Method of Solution We will use Matlab39s stiff solver called ode15s Also see the link to the ODE suite on the homepage for Math 302 Matlab Implementation The mfiles for this system are called yporegm and 0regm function yporegyporeg t y yporegl 7727y2y1y2y1 837510 6y1y1 yporeg2 y2 y1y2y37727 yporeg3 161y1 y3 yporeg yporeg1 yporeg2 yporeg339 your name your student number lesson number clear tf 1000 yo 4 11 4 t y ode15s39yporeg390 tfyo semilOgYty1Ity2ty3 title39your name your student number lesson number39 xlabel39time39 ylabel39solutions one two and three39 solutions one two and three 0 o o 0 your name your student number lesson number 100 150 200 250 300 time 350 400 Next the initial conditions will be chosen closer to the steady state solution so that man 4809 y20 11 and y30 480 your name your student number lesson number solutions one two and three I I I I I I I I I 0 100 200 300 400 500 600 700 800 900 1000 time Differential Equation Model Stefan39s law of cooling for a glowing hot mass ut cusur4 u4 For large values of u this will be a sli differential equation Method of Solution Euler 39s method for Newton law of cooling is uk ukh calm uk Solve for uk1 to get uk1 1 hcuk husurc Need 1 ch gt O Eulertrapezoid method where the right side is now evaluated at the next time step For the Newton cooling problem this uk ukh Cusur uk Solve for uk1 to get uk ukl hc husurcl hc Matlab Implementation of Eulertrapezoid function feuls feulstx feuls 69411 x273A4 your name your student number lesson number clear Euler Trapezoid Algorithm with Picard Solver eps 0001 maxm 10 u1 973 T 100 KK 40 h TKK tl 0 Begin Time Loop for k lzKK oldfeuls feulstkuk oldu uk holdfeuls tk1 tk h Begin Picard Loop for mlzmaxm newu uk h5oldfeuls feulstk1oldu if absnewuoldulteps break end oldu newu end uk1 newu end plot tu title 39your name your student number lesson number 39 xlabel 39 time 39 ylabel 39 temperature 39 Numerical Experiments The following system for two functions has known solutions y1t eXp4t and y2t eXpt y139 10004 y1 10000 yz4 Y239 2 Y139 Y2 39 3 24 Matlab Solvers 0de45 0de23s and Others One can use either Matlab s odedemo buiode or ypcstiffm and cstiffm les Instructor Inclass Matlab Calculations 0 Consider the Van der P01 1 describe derivation 0 Examine the code le vdp01m 1 point out the input execution and output sections 2 indicate the input variable 0 Examine the function f11e ypvdpolm 1 discuss the three meanings of the symbols y1 and y2 name function or array 2 go over the two differential equation lines and point out p o Execute the code f11e vdp01m 1 experiment with ode45 ode23s and ode15s 2 illustrate the graphs for charge and current 3 observe the number of time steps for u 10 50 and 100 Instructor Inclass Matlab Calculations Consider the population model Examine the code le rk4m 1 point out the input execution and output sections 2 indicate the input variables Examine the function le fpopm 1 indicate the maximum population and growth terms 2 evaluate the function which requires variables Execute the code le rk4m Explain the graphical output and the numerical output use the command whos Use the command hold on to generate several curves on the same graph Do several executions of rk4m with variable initial populations Consider the population model with variable growth term c and indicate how rk4m and fpopm need to be changed Differential Equation Model Newton39s law of motion then gives the differential equation mutmgku2 As time evolves the speed will approach u mg10m Method of Solution Euler39s Numerical Method for Approximating y39 gty and y0 given Yi1 Yi dtgidtYi Improved Euler39s Numerical Method for Approximating y39 gty and y0 given Yi1 Yi dtgidtYi Yi1 Yi dt2gidtYigi1dtYi1 Table Discretization Errors KK erreul errimerr 8 31152 01370 16 15347 00326 32 07571 00080 64 03761 00020 Matlab Implementation function fmass fmassty fmass 32 1absy your name your student number lesson number clear Y1 0 T 100 KK 100 h TKK t1 0 for k lzKK yk1 yk hfmasstkyk tk1 tk h yk1 yk 5hfmasstkyk fmasstk1yk1 end speed plottY title39your name your student number lesson number39 xlabel39time39 ylabel39speed39 yournameyourstudentnumbenlessonnumber 0 100 150 200 250 300 350 l 0 The second numerical experiment is with the resistive force equal to 1absu 001u2 yournameyourstudentnumbenlesson number speed Lesson 4 Population Taylor and RungeKutta Methods 41 Applied Problem Consider a single sh population whose size is given by xt The change in the size of the fish population over a given time interval is approximately equal to the product of the change in time At the difference in the birth rate and death rate b d and the quotpresentquot fish population xt This may be written as xt At xt m At b d xQ where t is between t and t At The simplest model is for b d being a constant A second model is the logistic model where b d cM x with c M constant and M is the maximum population size A third model may involve harvesting ht 0r stocking st at a given rate We want to be able to predict the fish population size as a function of these parameters In this lesson we will focus on the logistic model and in the lesson five we will account for harvesting or stocking of the lake In both these lessons we have assumed there are no predator fish that significantly alter the fish population we will study this in lesson seven 42 Differential Equation Model A differential equation can be formed from the above approximation by dividing by At and letting it go to zero The approximation will get better as At goes to zero because the population size in the interval will vary less as At goes to zero The above approximation becomes x b dx If b d are constant then the solution is xt x0 eb 39 d If b d cM x then the solution can be found in terms of a more complicated expression with an exponential function provided c and M are constants In more sophisticated models this is not the case and so we are forced to consider numerical solution to these models A variation of the logistic model that includes harvesting is x cM xx ht 43 Taylor Numerical Method We wish to have a higher order approximation scheme and so we apply the fourth order Taylor polynomial approximation to ut h and a t Let u have five continuous derivatives so that for some C between t and t h we have ut h ut um h unt 1122 umt h36 uma 11424 u5C h5120 Since u is a solution of the differential equation one can compute the higher order derivatives by using the chain rule for example 11a tut 1 futut ut fttut 1 fut11t ftutD If ftu is complicated then these calculations will become costly However for simpler cases it is very powerful Consider the Newton39s law of cooling where ftu cusur u Then the n3911 derivative of u has the form cquot391fu The above Taylor polynomial with variable h is u t h ut fut h c1 fut 1122 c2 fut h36 c3 fut 11424 c4 fCh5120 This suggests the following algorithm k1 U uk folk h cfuk 1122 c2 folk h36 c3fuk 11424 Since the error in Taylor39s theorem is of order 5 the discretization error might have an error of order 4 This is the case provided the solution has 5 continuous derivatives The following Matlab code in the mfile tayerrm illustrated the order of convergence for the cooling problem with c 126 ysur 70 and y0 200 The errors are computed at time equal to 100 and the steps sizes are reduced by 50 of the previous step size The table indicates the fourth order Taylor method has an error that is proportional the step size raised to the fourth Matlab m le tayerrm Euler and Fourth Order Taylor Errors clear kk 40 T 100 dt Tkk u0 200 c 126 usur 70 uk uO ut uO for k kk uk uk dtcusur uk uex usur u0 usur exp ckdt utp ut ut ut dtc usur utp ut ut th22ccusur utp ut ut th36cA2cusur utp ut ut th424cA3cusur utp erreul absuk uex errtay absut uex end kk erreul errtay Table Discretization Errors with Euler Order One and Taylor Order Four K h 10K Euler Error Taylor Error 5 2000 13696 10589 107 10 1000 06693 64089 107 20 0500 03309 03941 107 40 0250 01645 00244 107 The 11th order Talyor method for approximating the solution of ut ftut with u0 given is uk1 uk ftk uk h f0 tk uk hzZ fln391 tk uk hquotn where f6 is the ith derivative of the composed function ft ut Example In the above Newton39s model of cooling the derivatives of the composed function are f0 cifu provided fu cusur u with usur equal to a constant Example In the calculations below we consider the second order Taylor method but for more general ft ut uk uk ftk uk h m tk uk 1122 uk ftk uk h 110 uk fuck uk ftk uk 1122 44 RungeKutta Numerical Method The Taylor method for more accurate numerical solutions of ordinary differential equations requires the computation of a number of partial derivatives of ftu and this made its implementation difficult for higher order methods Here we will show how one can make use a number of function evaluations of just ftu and still obtain higher order methods based on the Taylor method The second order Taylor algorithm is uk uk h ftk uk tk uk fuck uk f tk ukh2 In order to avoid some of the partial derivatives we may evaluate ftu at a number of points The location of these points can be found by using the two variable version of Taylor polynomials for fXy The derivation can be done via the one variable problem where Ft fathbtk Use the extended mean value theorem F01 FO F1010 F2c1022 where c is in 01 The first and second order derivatives of Ft must be computed via the chain rule Fl thfykthDykf F Dx h Dy kzf Dx h Dy kfx h Dx h Dy kfyk fxxhh 2fxyhk fyykk Therefore the second order Taylor polynomial with error in two variables is fahbk fab fxab h fyab k fxxachbckhh 2fxyachbckhk fyyachbckkk2 This multilinear approximation of a function of two variables is the model which allows us to simplify the Taylor method for the numerical solution of ordinary differential equations Consider the above second order Taylor method We want to replace the expression in the square brackets with a suitable evaluation of f at ftk ac uk 5 The ac and 5 can be chosen via the second order Taylor polynomial error in two variables withatkh0cbukandk ftk ac uk a ftkuk tkuk x futkuk s O0c2oc z In order to match this with the terms in the square bracket in the second order Taylor algorithm let 0c h2 and l ftkuk h2 so that a variation of this method is uk uk h ftk m uk ftkuk h2 This is known as the midpoint method or the second order Runge Kutta method The fourth order RungeKutta method can be formulated using a similar construction It also can be viewed as Simpson s rule numerical approximation of the equivalent integral equation ut h ut T z uzdz t The integral can be approximated by h6k1 4kz2 k32 k4 where k1 ftut k2 ftl 2 ut h2k1 k3 fth2 ut h2kz and k4 fth ut 1113 The fourth order Runge Kutta method where t and ut are replaced by tk and uk is uk uk h6k1 2k2 2k3 k4 In order to show this scheme is fourth order convergent that is uth uk1 I S Ahs one must use the Taylor expansions and in this case the solution should have continuous derivatives of order ve 45 Matlab Implementation The following Matlab code for the RungeKutta fourth order method uses fixed step sizes The logistic differential equation is solved where c 01 M 100 and XO 10 Notice the solution approaches 100 The steady state solutions are defined by any constant where the right side is zero This requires cXM X 0 and so X 0 or X M The right side of the differential equation is stored is the function file fpopm and the RungeKutta code is stored in the file rk4m This initial population is stored in Xl and the values for c 001 and M 100 are in fpopm function fpop fpoptx fpop 01100 xx your name your student number lesson number clear x1 10 T 10 KK 50 h TKK t1 0 for k 1KK k1 hfpoptk xk k2 hfpoptk5hxk5k1 k3 hfpoptk5hxk5k2 k4 hfpoptkhxk5k3 xk1 xk k1 2k2 2k3 k46 tk1 tk h end plottx title your name your student number lesson number xlabel time ylabel population your name your student number lesson number 100 m o l 01 o I population 46 Numerical Experiments In the following numerical experiments we keep c 001 and M 100 fixed and vary the initial population We have executed rk4m with a number of initial conditions stored the output and plotted then all on the same graph Notice if xO is positive then the solution will always converge to the largest steady state solution Unless xO 0 the smallest steady state solution the xt always goes away from the smallest steady state solution Here the largest steady state solution is called stable and the smallest steady state solution is call unstable Use the hold on command and siX different initial populations 1 30 60 90 110 and 140 These should be stored in the array X ofthe rk4m le This is execute rk4m siX times with different Xl Xl l Xl 140 Note the abuse of notation if X Xt is a function oft then Xl is a function evaluation at t 1 if X is an array then Xl is the first entry in the array 140 120 O O I O populations for different x0 47 Additional Calculations Next we x the initial condition x0 10 M 1001S and vary the c You may want to increase the nal time from T 10 to T 20 or larger a State the logistic differential equation for c 0005 and M 1001S b Modify the popm and rk4m les c Execute the rk4m le for c 0005 d Use the hold on command and repeat c for c 001 0015 and 002 e Examine the curves What happens to the population as c increases Consider shhm and yp shhm Fix the initial condition X0 10 M 100lS ht 401sinpit and vary the c You may want to increase the nal time from T 10 to T 20 or larger a State the logistic differential equation for c 0005 and M 100lS b Modify the ypflshhm and f1shhm les c Execute the f1shhm le for c 0005 d Use the hold on command and repeat c for c 001 0015 and 002 Examine the curves What happens to the population as c increases Differential Equation Model Let S 7727 W 161 q 837510A6 y10 4y20 11 and y30 4 3 139 sY2 39 YI Y2 YI 39 q 3 12 Y239 y2 y1 y2 y3s and y339 WYl 39 Y3 Matlab39s 0de23s and 0de15s can be used to solve such systems The steady state solutions which are defined by all three derivatives set equal to zero so that 0 sY2 39 YI Y2 YI 39 q 3 12 f1 h Y2 Y3 a 0 39 Y2 39 YI Y2 Y35 f2 h Y2 Y3 and 0 WYI 39 Y3 f3 1a Y2 Y3 The three positive steady state solutions are ylquot y 48818 and y 99796 These may be found by using the symbolic Matlab command solve as follows EDU y1y2solve39y2 y1y2y1 qy1y103939y2 y1y2y1039 The Jacobian of f1 f2 f339 is f1 f1 fa sy2 1q2y1 51y1 0 J f2M f2 f2 y2s 1 y1s 1s f3M f3h Eh w 0 w The eigenvalues of J evaluated at the steady state solutions are easily found by using the Matlab command eig The following code is in the Matlab m le called 0regjacm yl 48818 y2 99796 Y3 yl s 7727 w 161 q 8375lOA 6 J 8 y2l Q2yl 81 yl 0 y2s l yls 15 W 0 w eigJ oregjac ans 257146 187473 00013 Method of Solution We will use Matlab39s stiff solver called ode15s Also see the link to the ODE suite on the homepage for Math 302 Matlab Implementation The mfiles for this system are called yporegm and 0regm function yporegyporegty yporeg1 7727y2y1y2y1 837510 6y1y1 yporeg2 y2 y1y2y37727 yporeg3 161y1 y3 yporeg yporeg1 yporeg2 yporeg339 your name your student number lesson number clear tf 1000 yo 4 11 4 t y ode15s39yporeg390 tfyo semilogYtylIty2ty3 title39your name your student number lesson number39 xlabel39time39 ylabel39solutions one two and three39 solutions one two and three your name your student number lesson number 0 50 n I n n 100 150 200 250 300 350 400 Next the initial conditions will be chosen closer to the steady state solution YI0 4809 y20 11 and y30 480 your name your student number lesson number 10 solutions one two and three 0 100 200 300 400 500 600 700 800 900 1000 time Instructor Inclass Matlab Calculations 0 Consider the massspring problem 0 Describe the subroutine ode45 for systems 0 Examine the code le msm 1 point out the input execution and output sections 2 indicate the input variables 3 indicate the initial array 0 Examine the function f11e ypmsm 1 discuss the three meanings of the symbols y1 and y2 name function or array 2 go over the two differential equation lines o Execute the code f11e msm 1 explain the graphical output 2 discuss the initial array and use variable initial conditions 3 use variable c 005 01 and 05 4 use variable w 08 10 and 12 in ft sinwt Differential Equation Model View the string as three masses lumped at the intersections of segments and with displacements uit where i l 2 and 3 The mass of each segment is p AX so that mass times acceleration is p AX ui There are three forces external from the pressure AX ft from the tension on the right Tui1 ui AX and from the left Tui ui1 AX From Newton s law of motion we have p AX ui AX ft Tui1 ui AX Tui ui1 AX Divide by p AX and define OL Tp AX2 to get the system of three second order differential equations ODE Model for Discrete String 11in f OLLli1 L1H 20L Lli ui0 O and ui 0 O 2 u0t unt O for t gt O 3 Equation 1 may be put into the matrix version of a system of ODEs For example if the string is divided into four equal parts then 11 4 and 1 may be written as siX scalar equations 111 2 114 L12 2 L15 L13 2 L16 u4 u1 lpcf OLL12O 20111 1153 2 L12 2 f OCL11 L13 20C L12 1163 2 L13 2 f OLO L12 20C L13 Or using a vector equation and a 6X6 matrix in u ul sin96812t u u2 J5 sin96812 r u39 3M u3 abzap sin96812t and u4 M4 0 M5 0 L16 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 2 Aza 2 1 0 0 0 0 w1thaTpAx 1 2 1 0 0 0 0 1 2 0 0 0 Method of Solution Matlab s ode45 can be used to nd numerical solution of system of siX ODEs The numerical method is a variation of the RungeKutta 45 variable time step method Matlab Implementation Since there are three unknowns and three second order differential equations and we wish to use Matlab39s ode45 scheme the file ypstringodem must contain the three right sides of the differential equations where y1 displacement of the left segment y2 displacement of the center segment and y3 displacement of the right segment y4 velocity of the left segment y5 velocity of the center segment and y6 velocity of the right segment The stringodem file contains the call to ode45 and the Matlab command plot generates the graphs of the three displacements The initial displacements are 10 20 20 and initial velocities are 00 00 00 They are stored in the uo array in the stringodem le mfiles ypstringodem function ypstringode ypstringodetu T 10 rho 1 L 1 dX L4 cc TrhodXdX 10sin96812t1 1414 10 O Orho ypstringode1 u4 ypstringode2 u5 ypstringode3 u6 ypstringode4 cc2u1 u2 f1 ypstringode5 ccu1 2u2 u3 f2 ypstringode6 cc2u3 u2 f3 ypstringode ypstringode1 ypstringode2 ypstringode3 ypstringode4 ypstringode5 ypstringode639 mfiles stringodem clear clf T 10 rho 1 L 1 dX L4 cc TrhodXdX uou2200m to O tf 8 t u ode4539ypstringode39to tfuo man sizet1 X O251 gure1 for k 11man plotXO uk1 3 0 axis0 1 6 6 tit1e39string at various times39 Xlabe139position on string39 hold on pause end Am001oq OOOOIQ oooooh 2cc cc 0 O O 0 cc 2cc cc 0 O O 0 cc 2cc O O O eigenvalues eigA gure2 plottu1tu2tu3 Xlabe139time39 The rst graphical output has multiple curves each associated with a particular time times eigenvalues 0 233725 03098 178885 10076 96812 17680 233725 25689 178885 33798 96812 42189 50934 59623 68554 77659 The amplitudes increase as time increases This is also depicted in the second graph which has the displacements of the three segments versus time The increasing amplitudes occur because the forcing term is f10sin96812t1 1414 1 0 0 0rho and has the same frequency as associated with the smallest positive eigenvalue The third gure uses 86912 in f not 96812 and indicates the discrete string is not resonating string at various times i i i i i 0 01 02 03 04 05 06 07 08 09 1 position on string x x x w I w I w I x 8 0 1 2 3 4 5 6 7 8 3 0 1 2 3 4 5 6 7 8 Differential Equation Model Stefan39s law of cooling for a glowing hot mass ut cusur4 u4 For large values of u this will be a sli differential equation Method of Solution Euler39s method for Newton law of cooling is uk1 ukh Cusur uk Solve for uk1 to get uk1 1 hcuk husurc Need 1 ch gt O Eulertrapezoid method the right side is now evaluated at the next time step For the Newton cooling problem this uk1 ukh cusur uk kl Solve for u to get uk1 ukl hc husurcl hc For the Stefan cooling problem we get a fourth order polynomial equation to solve a each time step ukl uk Cusur4 ukl4 Matlab Implementation of Eulertrapezoid function feuls feulstx feuls 69411 x273A4 your name your student number lesson number clear Euler Trapezoid Algorithm with Picard Solver eps 0001 maxm 10 u1 973 T 100 KK 40 h TKK tl 0 Begin Time Loop for k lzKK oldfeuls feulstkuk oldu uk holdfeuls tk1 tk h Begin Picard Loop for m1maxm newu uk h5oldfeuls feulstk1oldu if absnewuoldulteps break end oldu newu end uk1 newu end plottu title39your name your student number lesson number39 xlabel39time39 ylabel39temperature39 Numerical Experiments The following system for two functions has known solutions y1t eXp4t and y2t eXpt y139 10004 y1 10000 y24 3 239 2 Y1 39 3 2 39 3 24 Matlab Solvers 0de45 0de23s and 0de15s One can use either Matlab s odedemo buiode or ypcstiffm and cstiffm les Lesson 15 Oregonator Chemical Reactions and ode15s 151 Applied Problem The chemical reaction with bromous acid bromide ion and cerium ion exhibits a remarkable chemical attributes which oscillate with changes in color and structure The rate of changes of the concentrations are very large and have spikes similar to those in the Van der Pol equation solution Such systems are called stiff and require special numerical methods The Matlab command odedemo illustrates this see chm9ode and other stiff differential equations and two stiff solvers called ode23s and ode15s 152 Differential Equation Model The chemical reaction can be described in five stages with known reaction rates Upon suitable time and concentration scaling this is modeled by three coupled differential equations Let s 7727 w 161 q 837510A6 y10 4y20 11 and 3130 4 W sy2 y1y2 yl q ylz Y239 Y2 yl Y2 f Y3S and y339 Wyl Y3 Matlab s ode23s and odelSs can be used to solve such systems As a first step in the study of this system we consider the steady state solutions which are defined by all three derivatives being set equal to zero so that 0 SY2 y1y2 y1 q y12 f1y1 yz Y3 0 Y2 yl Y2 f Y3S f23 1 Y2 Y3 and 0 WWI Y3 fsY1 Y2 Y3 The three positive steady state solutions are ylquot y 48818 and yzquot 99796 These may be found by using the symbolic Matlab command solve as follows EDU y1y2solve39y2 y1y2y1 qy1 y103939 y2 y1 y2y1039 y1 W 12qqqA28IIA12l 12qIIqA28IIA12l y2 0l 114q 14qA28qA12 114q14q 28q 12 The Jacobian of f1 f2 f3 is fly flyz fly Sy2 1 qzyl 51 yl 0 J flyi fin fly 2 yZS 1Y1S 15 3m 3y2 f3 W 0 W The eigenvalues of J evaluated at the steady state solutions are easily found by using the Matlab command eig The following code is in the Matlab m le called 0regjacm y1 48818 y2 99796 Y3 y1 s 7727 w 161 q 837510A 6 J s y21 q2y1 s1 y1 0 y2s 1y1s 15 w 0 w eigJ The ouput from the execution of this le shows some large positive and negative eigenvalues so that the steady state solutions are unstable and the system is stiff near the steady state solutions EDU oregjac ans 257146 187473 00013 153 Method of Solution We will use Matlab39s stiff solver called ode15s This is much more sophisticated than the simple Eulertrapezoid method that was used a previous lesson However its formulation is somewhat similar to the implicit nature of each time step in the Euler trapezoid method Further details about this method are described in the Matlab s helpdesk also see the link to the ODE suite on the homepage for Math 302 154 Matlab Implementation The mfiles for this system are called yporegm and oregm In the first calculation we used the final time tf 300 and the second calculation used tf 1000 so as to show the oscillations over three periods The graph is for the scaled concentrations versus time The oregm file uses semilogy graphing and so it reveals very signi cant spikes in the solutions function yporegyporegty Yporeg1 7727Y2Y1Y2Y1 83751OA6Y1y1 Yporeg2 Y2 Y1Y2y37727 yporeg3 161y1 y3 yporeg yporegm yporeg2 yporeg3139 your name your student number lesson number clear tf 1000 yo 4 11 4 t y ode15s yporeg 0 tfyo semilogytYI1tY2tY3 title your name your student number lesson number xlabel time ylabel solutions one two and three your name your student number lesson number solutions one two and three your name your student number lesson number two and three 3N 0 solutions one I I I I I I I I I 0 100 200 300 400 500 600 700 800 900 1000 time 155 Numerical Experiments In the calculation below we decreased w from 161 to 04 Notice the period increased and the steepness of the spikes decreased your name your student number lesson number o o solutions one two and three I I I I I I I I I 0 100 200 300 400 500 600 700 800 900 1000 time Next the initial conditions will be chosen closer to the steady state solution so that y10 480 y20 11 and y30 480 The rst calculation reveals the same oscillating pattern is established The second calculation is for a much smaller time interval where tf has been decreased from tf 1000 to 10 your name your student number lesson number two and three solutions one I I I I I I I I I 0 100 200 300 400 500 600 700 800 900 1000 time your name your student number lesson number 10 l I I I I l I l I 0 solutions one two and three 156 Additional Calculations Consider the above oregonator chemical reaction with variable w S10 S and 2S a State the new system of differential equations b Modify the yporegm le c Insert your name and student number into the oregm le Execute the oregm le for the three choices of w d Print the yporegm and oregm les and the graphs A Lesson 4 LR Circuit and Variable Time Steps 41 Applied Problem In the previous sections we introduced the improved Euler method to obtain more accurate numerical solutions of ordinary di erential equations m ft ut and u0 given These schemes required the user to select a xed step size h If the solution has a large derivative then smaller step sizes will be needed to track the solution If the solution is almost linear then much larger steps can be used Since many problems have both types of solution segments it is desirable to change step size while the scheme is being executed Here we will show how second order methods can be used to adjust the step size of the first order Euler method We will consider a series LR circuit where the imposed voltage has a high frequency and therefore we may require variable time steps In a later lesson we will use a much more sophisticated method which also uses variable step sizes to study tuning a series LRC circuit 42 Differential Equation Model Consider a di erential equation whose right side signi cantly varies with time This means the derivative of the solution must have large changes with time One example is the Newton cooling model where the surrounding temperature is not constant A related example is a series LR circuit with time dependent imposed voltage LI RI v v0 Sinot Here I is the current L is the inductance R is the resistance and V is the imposed voltage This differential equation is derived from Kirchhoff39s law which requires the voltage drops over each component to sum to the imposed voltage The voltage drop over the inductor is LI and the voltage drop over the resistor is RI 43 Method of Solution The local error in a numerical method is uk1 utkh where h is the step size In rst order methods this is bounded by Ahz and in second order methods this is bounded by Bh3 The constants A and B will depend on some bounds of higher order derivatives Reasonable estimates for A and B can be computed if the exact solution utkh is known but most of the time that is what we are trying to nd The exact solution will be modeled by the higher order method Uk where IUkH utkh s Bh3 1 The lower order method uk satis es luk utkh s Ahz 2 In our calculations we will know both uk1 and Uk but not utkh Now write uk Uk uk utkh utkh Uk Then by the triangle inequality and the error bounds l and 2 we have 11k Uk s 11k utkh utkh Uk 3 A112 13113 sAh21BAh BA h is the error made by replacing utkh with U Thus we may approximate A by A luk1 Uk1 lhz Furthermore the local error for the lower order scheme can be approximated by luk utkh In Ahz 11k Uk Ihz hz 3 Let tol tolerance be a positive real number for which we want the error to be less than Then from 3 we should choose h so that A 112 3 m1 At each attempted time step we check to see if this is close to being satis ed If A h2 is much less than tol then we may consider making h larger Or if A h2 is much larger than tol then we should consider making h smaller The new value of h is o en viewed as a scalar multiple of the old h that is newh scaleh The scaling factor can be found from A newh2 tol A scalequot h2 tol Then scale tol I uk1 Uk1 I 12 44 Matlab Implementation Another variable step method is to formulate his in terms of the local truncation error terr uk utkhh This is done in the code below where scale tolterr The following is a slight variation on the Matlab code for Euler s method There is a single loop which corresponds to time steps This code approximates the solution to u 570 30sin11t u and u0 200 The second order improved Euler method is used to estimate the exact solution and then used to nd the new step size for the rst order Euler method Euler newy Y hftY Improved Euler newY Y h2 ftY fthnewy Approximate Truncation Error newY newyh Y h2 ftY fthnewy Y hftY h2 ftY fthnewy In the code below we consider the cooling problem where the surrounding temperature varies by a trig function The right side of the di erential equation is stored in the m le called rm function fl fl t y flr 2 70 30sin pit 39 Y The variable step size code using Euler and improved Euler methods is stored in the m le called eulervm Notice the for loop command is replaced by the while loop command This is necessary because we desire the simulation to run over a given time interval here from time equal 0 to 20 but the time step sizes will be changing according to the scale computation at each time step Note the if statement requires terrlttol 11 and the 11 factor is inserted so that the following statements for adjusting the h do not create an in nite loop In order to observe such and undesirable event replace the 11 by 10 and use tol 10 your name your student number lesson number clear y1 200 tol 5 h 01 t1 0 k 1 while tklt20klt2000 newt tk h oldfff flrtkyk newy yk holdfff terr abs5 oldfff flrnewtnewy if terrlttol11 tk1 newt yk1 newy k k1 end scale tolterr if scalelt1 h 1h end if scalegt10 h 10h end if scalegt1scalelt10 h scaleh end end steps sizet2 dt difft plott1steps 1y1steps 1t1steps 1300dt title your name your student number lesson number xlabel time ylabel current and time step size EDU eulerv steps 270 In the graphical output below the top curve is the approximate solution given by Euler39s method with variable time steps The lower curve is an ampli ed version of the change in time step versus time Note the time steps are much larger when the solution is almost linear There are a number of potential problems If tol is too small then the h will be very small and this will result in too many time steps If tol is too large then h will be large and the local errors may not be acceptable Ifthe calculation of A gives a real small number then there may be signi cant roundolT errors Generally one should put lower and upper bounds on the computed scaling factors and step sizes your name your student number lesson number current and time step size 45 Numerical Experiments The tol parameter in the above calculations was equal to 5 and it took 270 time steps to reach the nal time equal to 20 Ifthe tol parameter is set equal to 10 then it will require 150 time steps For the tol parameter equal to 1 it takes 1134 time steps The above code is easily broken Fortunately there are a number more robust variable step codes for example see the Matlab applications ode23 and ode45 Consider the LR circuit with 10 0 L 12 R 10 and V 30 sin1110t In eulervm use tol 5 and use the nal time equal to 2 in the while command your name your student number lesson number current and time step size 46 Additional Calculations Consider the LR circuit with 10 0 L 12 R 2 and V 30 sin1130t In eulerVm use tol 10 and use the nal time equal to 4 in the while command a State the diiTerential equation b Modify the rm 0 Insert your name and student number into the eulervm le and modify this le for the new tol and nal time Execute the code in the eulervm le d Print the rm and eulervm as well as the graph Differential Equation Model Series LR circuit with time dependent imposed voltage is LI39 R1 V V0 sinoat L inductance 1 current R resistance and V imposed voltage Method of Solution Use variable step sizes We would like to adjust the step size dt h so that error uk1 uklh is small In general we do NOT know uklh this is the reason for using numerical methods Therefore we do not know the error One way around this is to replace the exact solution by a higher order numerical solution For example in order to estimate the error for Euler s method approximate the exact solution by the Euler trapezoid method l Uk1 utkh l 3 13h3 l uk1 utkh l g Ahz In our calculations we will know both uk1 and U but not utkh uk1 Uk1 uk1 utkh utkh Ukr uk Uk1 S uk utkh utkh Uk1 3 Ah2 Bh3 gAh2 1 BA h BA h is the error made by replacing utkh with Ukll Thus we may approximate A by A uk Uk1 lhz Furthermore the local error for the lower order scheme can be approximated by ukl utkh z A 112 ukl Ukl hZ Let tol tolerance be a positive real number for which we want the error to be less than We should choose h so that A h2 S tol Matlab Implementation Euler newy Y hft7Y Improved Euler newY Y 112 ftY fthnewy Approximate Truncation Error newY newyh Y h2 ftY fthnewy Y hftY 112 ftY fthnewy Use the function le rm and the m le eulervm function flr flrty flr 270 30sinpit y your name your student number lesson number clear Y1 200 tol 5 h 01 tl 0 k 1 while tklt20klt2000 newt tk h oldfff flrtkyk newy yk holdfff terr abs5oldfff flrnewtnewy if terrlttol11 tk1 newt yk1 newy k k1 end scale tolterr if scalelt1 h 1h end if scalegt10 h 10h end if scalegt1scalelt10 h scaleh end end steps sizet2 dt difft plott1steps1y1steps lt1stepsl300dt title39your name your student number lesson number39 xlabel39time39 ylabel39current and time step size39 current and time step size your name your student number lesson number If tol gets small then the time step will get smaller and more steps will be required to get to the nal time For example Tol 10 requires 150 time steps Tol 5 requires 270 time steps and T01 1 requires 1134 time steps Homework and Report for Lesson 15 Your Name Your Class Number Section Number Lesson Number 1 Describe the Applied Problem 1 point 2 State the Differential Equation Model 1 point 3 Describe the Numerical Method 1 point 4 Include the Matlab Code 2 points 5 Numerical Experiments Consider the Oregonator chemical reaction system with variable w In your calculations be sure to use a suitably long time interval so as to see the periodic nature of the solutions a 3 points Find the graphical solution when W S10 b 3 points Find the graphical solution when w2S c 3 points Find the graphical solution when w 3S d 3 points What happens to the solution when w increases e Extra Credit for 3 points Modify the light side of the third equation to y339 wy1 y3 1 Use ode15s to find the solution What are the steady state solutions Are they stable Consider eulervm and rm Consider the LR circuit with 10 0 L 12 R 2 and V 30 cos730t In eulerVm use tol 10 and use the nal time equal to 4 in the while loop a State the differential equation b Modify the rm c Modify eulerVm for the new tol and final time Execute the code in the eulerVm file Use the whos command to record the number of time steps used d 6 Repeat c for tol 5 and tol 1 Compare the graphs generated using tol 10 5 and 1 Examine the graphs What happens as tol decreases
Are you sure you want to buy this material for
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'