## Transcribed Text

Computational Math
For this programming assignment you may use Matlab. Your tasks are:
(4.1.a) to implement the 6 stage Fehlberg embedded Runge Kutta 4,5 pair
(4.1.b) to demonstrate their correctness integrating problems that should produce “exact” results and others with known solutions;
(4.1.c) to investigate the performance of the order 4 and order 5 methods run individually with fixed time steps versus the performance of the pair together with
automatic stepsize control. (You should use equations of your choice but must
include at least the example presented below.)
The class webpage contains a Matlab code that uses an Adams Bashforth-Moulton integrator (ode113) from Matlab to integrate an astronomy example from Asher and Petzold’s
text book. Your explicit RK embedded pair code will not produce exactly the same trajectory due to the singularity of the differential equation around (ˆµ, 0) unless you modify the
code that evaluated the derivatives, however, it should look close to that produced by the
example. The system of differential equations integrated in this example is
D1 = ((y1 + µ)
2 + y2
2)
3/2
D2 = ((y1 − µˆ)
2 + y2
2)
3/2
y
1 = y3
y
2 = y4
y
3 = y1 + 2y4 − (ˆµ((y1 + µ)/D1)) − (µ((y1 − µˆ)/D2))
y
4 = y2 − 2y3 − (ˆµ(y2/D1)) − (µ(y2/D2))
y1(0) = 0.994
y2(0) = 0.0
y3(0) = 0.0
y4(0) = −2.00158510637908252240537862224
µ = 0.012277471
µˆ = 1 − µ
The equation should exhibit periodic behavior and integrating on the interval [0, 17.1] should
produce the trajectory of interest in the y1, y2 plane.
When integrating with a fixed stepsize explicit RK method it may take many steps. You
should set the stepsize to use 100 steps, 1000 steps, 10000 steps etc. until you see behavior like
1
that produced by the example code or until you are sure it cannot be reproduced efficiently
using a fixed stepsize with either the order 4 method or the order 5 method in your code.
You should then run your automatic stepsize code to see if you can produce the desired
trajectory (qualitatively at least) in a number of steps comparable to the Matlab integrator
ode45 which is also an embedded RK pair with automatic stepsize selection.
Stepsize Selection: An embedded pair gives yn and ˜yn that are order p and order p + 1
respectively. Therefore, n ≤ Chp+1 ≈ |y˜n − yn| for yn. Given error tolerances we have
|y˜n − yn| ≤ τrel|yn| + τabs → accept the step, otherwise reject the step and reduce stepsize
|y˜n − yn| ≤ µ1(τrel|yn| + τabs), 0 < µ1 < 1 → accept the step and consider increasing the stepsize
Note that these conditions can be imposed on each component of y individually or using a
vector norm.
The stepsize reduction finds α < 1 such that
αp+1|y˜n − yn| ≤ µ2τrel|yn| + τabs, µ1 ≤ µ2 < 1
where µ2 is a safety factor to make it unlikely the new stepsize will be rejected. When
increasing the stepsize 1 <β< 1 + µ3 such that
βp+1|y˜n − yn| ≤ µ2τrel|yn| + τabs,
so h cannot increase too quickly.

These solutions may offer step-by-step problem-solving explanations or good writing examples that include modern styles of formatting and construction
of bibliographies out of text citations and references. Students may use these solutions for personal skill-building and practice.
Unethical use is strictly forbidden.

% Usage: [y tme] = rkf45_mod(f,a,b,ya,h,rtol,atol,flg)

% Runge-Kutmea-Fehlberg 4th Order/5th Order embedded pair

%

% Input:

% f - ODE function f(t,y)

% a,b - Time interval of integration

% ya - initial condition

% h - initial step-size

% rtol - relative error tolerance per step

% atol - Absolute Error Vector

%

% Output:

% y - computed solution

% tme - time vector

%

% Examples:

% [y t]=rkf45(@myode,0,10,[1 0 0],1E-2,1e-1,1e-4); -Here, 'myode' is a user-defined function in M-file of the required ODE system

function [y,tme] = rkf45_mod(f,a,b,ya,h,rtol,atol,flg)

%Fehlberg RK45 Butcher Tableau coefficients

c10 = 0;

c20 = 1/4; c21 = 1/4;

c30 = 3/8; c31 = 3/32; c32 = 9/32;

c40 = 12/13; c41 = 1932/2197; c42 = -7200/2197; c43 = 7296/2197;

c50 = 1; c51 = 439/216; c52 = -8; c53 = 3680/513; c54 = -845/4104;

c60 = 1/2; c61 = -8/27; c62 = 2; c63 = -3544/2565; c64 = 1859/4104; c65 = -11/40;

%4th order coefficients

cw1 = 25/216; cw2 = 0; cw3 = 1408/2565; cw4 = 2197/4104; cw5 = -1/5; cw6 = 0;

%5th order coefficients

cz1 = 16/135; cz2 = 0; cz3 = 6656/12825; cz4 = 28561/56430; cz5 = -9/50; cz6 = 2/55;

%Error Coefficients

ce1 = 1/360; ce2 = 0; ce3 = -128/4275; ce4 = -2197/75240; ce5 = 1/50; ce6 = 2/55;

% Numerical Parameters

alph = 0.8; k = 0;

% Initial time moment

i = 1; tme(1) = a; t = a;

% Initial condition

y(1,:) = ya; y0 = ya;

% If it is the last iteration, then l_it = 1, othery0se l_it = 0

l_it = 0;

if flg==1%Perform 4th Order RK

while l_it == 0...