Complete the following problems with a single script and custom functions as appropriate. Important lines of code should include descriptive comments.
AERO 300 Laboratory Runge-Kutta-Fehlberg Methods The pre-lab assignment (Section 4) is due at the beginning of this lab (can be hand-written or typed). The lab (Section 5) will be due before your next lab section (.zip file submitted to PolyLearn). 1 Objectives This lab is intended to provide an overview of higher order Runge-Kutta methods for solving ordinary differential equations. The topics covered are: • Runge-Kutta-Fehlberg method • ode45 At the completion of this lab you should be able to solve ordinary differential equations using these tools. 2 Introduction By now, you should be familiar with the Runge-Kutta method for solving ordinary differential equations and initial value problems. As discussed in class, there is a local and global error associated with these methods which is related to the user defined step size, ℎ. A smaller value of ℎ will result in a more accurate solution at the expense of additional function evaluations. Therefore, a reduced step size results in a longer computational time. E. Fehlberg proposed and developed an error control method by using two Runge-Kutta methods of different orders to go from (?? , ??) to (??+1, ??+1). The difference of the computed y-values at ??+1 results in an estimate of the local error that can be used to control the step size. This concept is known as an embedded pair. 2.1 Runge-Kutta Methods In lecture, we reviewed the Euler Method and the Explicit Trapezoid method. Both methods fall under the umbrella of Runge-Kutta methods. In general, the second-order Runge-Kutta methods can be expressed as in equation (6.49) from Sauer which is repeated below: ??+1 = ?? + ℎ (1 − 1 2? ) ?(?? , ??) + ℎ 2? ?(?? + ?ℎ, ?? + ?ℎ?(?? , ??)) We can derive higher order (order ?) Runge-Kutta methods using the following: ??+1 = ?? + ℎ ∑ ???? ? ?=1 + ?(?ℎ+1) where ?? = ? (?? + ??ℎ, ?? + ℎ ∑ ????? ? ?=1 ) The values of ??, and ???, can be chosen depending on where you want to evaluate the derivative on the interval. For more, see the Wikipedia article. 2.2 Runge-Kutta-Fehlberg Method Fehlberg discovered two Runge-Kutta formulas that together need only 6 function evaluations per step and provide a 5th order estimate of the error for adaptive step size control. These equations are also in Sauer but listed here for convenience (equations 6.62, 6.63, and 6.64). ?1 = ?(??, ??) ?2 = ? (?? + 1 4 ℎ, ?? + 1 4 ℎ?1) ?3 = ? (?? + 3 8 ℎ, ?? + 3 32 ℎ?1 + 9 32 ℎ?2) ?4 = ? (?? + 12 13 ℎ, ?? + 1932 2197 ℎ?1 − 7200 2197 ℎ?2 + 7296 2197 ℎ?3) ?5 = ? (?? + ℎ, ?? + 439 216 ℎ?1 − 8ℎ?2 + 3680 513 ℎ?3 − 845 4104 ℎ?4) ?6 = ? (?? + 1 2 ℎ, ?? − 8 27 ℎ?1 + 2ℎ?2 − 3544 2565 ℎ?3 + 1859 4104 ℎ?4 − 11 40 ℎ?5) ??+1 = ?? + ℎ ( 25 216 ?1 + 1408 2565 ?3 + 2197 4104 ?4 − 1 5 ?5) With embedded pair, ??+1 = ?? + ℎ ( 16 135 ?1 + 6656 12825 ?3 + 28561 56430 ?4 − 9 50 ?5 + 2 55 ?6) Therefore, the error for step size control is, ??+1 = |??+1 − ??+1| = ℎ | 1 360 ?1 − 128 4275 ?3 − 2197 75240 ?4 + 1 50 ?5 + 2 55 ?6| If the relative tolerance is defined by the user as, ?, then the relative error test is given as, ??+1 |??+1| < if="" the="" test="" passes,="" then="" the="" (approximate)="" solution="" at="" +1="" is="" found="" and="" the="" solver="" moves="" to="" find="" the="" (approximate)="" solution="" at="" the="" next="" time="" step,="" +2.="" furthermore,="" the="" approximate="" solution="" +1="" at="" +1="" is="" replaced="" with="" the="" locally="" extrapolated="" version="" +1.="" if="" the="" test="" fails,="" the="" step="" size="" is="" reduced="" as="" in="" equation="" 6.57="" from="" sauer="" (and="" repeated="" and="" corrected="" below),="" https://en.wikipedia.org/wiki/runge%e2%80%93kutta_methods="" ℎ∗="0.8" (="" |??+1|="" )="" 1="" 5="" ℎ?.="" with="" this="" new="" time="" step,="" the="" method="" is="" applied="" again="" to="" verify="" that="" the="" computed="" error="" satisfies="" +1="" |??+1|="">< .="" if="" so,="" we="" take="" +1="??" +="" ℎ∗.="" however,="" if="" using="" the="" time="" step="" ℎ∗="" still="" yields="" +1="" |??+1|=""> ?, then one should reduce ℎ∗ by, for instance, repeatedly dividing it by 2 (and applying the scheme again) until the error condition is satisfied. Suppose this occurs for ℎ = ℎ̅ . Then, we take ??+1 = ?? + ℎ̅ . In this type of scheme, iteration are now performed within a while loop since the number of time steps is not known a priori. Only the integration interval is known a priori. 2.3 MATLAB and ode45 The MATLAB function ode45 is used to solve ordinary differential equations. ode45 is based on an explicit Runge-Kutta formula known as the Dormand-Prince pair (see example 6.22 in Sauer). In general, ode45 is the best function to apply as a first try for most initial value problems. 3 Help Potentially useful terms for this lab: • ode45 • for • function • while • exp • disp • plot3 • axis • grid • hold • legend • title • xlabel • ylabel 4 Pre-lab Assignment NOTE: This must be done before lab. If it is not completed, you will not be allowed into lab. Review the help file for the ode45 function. Write a Matlab script to solve the IVP ?′′ + ??′ + ? = 0 with ?(0) = 1, ?′(0) = 0. Submit only a published version of your code and a plot of the solution of the ODE. Andy G Cristales 5 Lab Assignment Complete the following problems with a single script and custom functions as appropriate. Important lines of code should include descriptive comments. You are encouraged to work in groups, however the code you submit must be your own. All graphs/figures/sketches should include a grid, labels, legend (if necessary), and the appropriate fontsize/linewidth/markersize. 1. Write a well-documented function to perform the Runge-Kutta-Fehlberg method on a first-order ordinary differential equation. It should have the form: [t,y] = rkf45(fun, tspan, y0, h, rTol) where fun is the function, tspan is the time span, y0 is initial condition, h is the initial step size, and relTol is the error tolerance. Consider the problem: ?′ = (? − ? − 1)2 + 2, ?(0) = 1 a. Solve for ?(?) using your RKF function with ℎ = 0.01 on the timespan ? = [0 ? 3 ] and relative tolerance, 10−6. Plot your results. b. Solve for y, using the MATLAB ode45 function with an absolute tolerance of 10−6. Plot your answer. c. Compare your results of (a) and (b) with the exact solution ?(?) = 1 + ? + tan ?. Plot the errors. 2. Use your rkf45 and the ode45 functions to solve for the solution to the chaotic system known as the Lorenz equations on the interval ? = [0 50]. You must use a function file to define the Lorenz Equations; no anonymous functions for this problem. Plot your results on a three-dimensional graph using the plot3 command. You will need to pass in a vector of functions and initial conditions. Choose initial conditions and a relative tolerance which are appropriate for an accurate solution. 6 Material to be Submitted Submit the code for the lab assignment (titled: yourname_Lab7) into your lab sections PolyLearn Assignment before the start of the next laboratory. Zip all your code into one file and put your name on every piece of code you turn in. Do not forget the master file. https://en.wikipedia.org/wiki/Lorenz_system#/media/File:A_Trajectory_Through_Phase_Space_in_a_Lorenz_Attractor.gif https://en.wikipedia.org/wiki/Lorenz_system#/media/File:A_Trajectory_Through_Phase_Space_in_a_Lorenz_Attractor.gif https://en.wikipedia.org/wiki/Lorenz_system 1 Objectives 2 Introduction 2.1 Runge-Kutta Methods 2.2 Runge-Kutta-Fehlberg Method 2.3 MATLAB and ode45 3 Help 4 Pre-lab Assignment 5 Lab Assignment 6 Material to be Submitted AERO300: Aerospace Engineering Analysis Lab 7 Grading Rubric Instructor will average each lab section scores to a common mean. This means it is ok if different TAs grade slightly differently, it’ll all be averaged out at the end. Problem 1 – 70 points code 40 a 10 b 10 c 10 Problem 2 – 30 points code 20 plot 10 Total 100 AERO300: Aerospace Engineering Analysis Lab 7 Grading Rubric