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...

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 e
or
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 e
or 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 e
or 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 um
ella of Runge-Kutta methods. In general, the second-order Runge-Kutta
methods can be expressed as in equation XXXXXXXXXXfrom 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 e
or 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 e
or 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 e
or 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 co
ected 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 e
or 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 e
or 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 e
or tolerance.
Consider the problem:
?′ = (? − ? − XXXXXXXXXX, ?(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.
. 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 e
ors.
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 Ru
ic

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 Ru
ic