QuestionQuestion

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

Solution PreviewSolution Preview

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

By purchasing this solution you'll be able to access the following files:
Solution.zip.

$100.00
for this solution

PayPal, G Pay, ApplePay, Amazon Pay, and all major credit cards accepted.

Find A Tutor

View available Administration - Other Tutors

Get College Homework Help.

Are you sure you don't want to upload any files?

Fast tutor response requires as much info as possible.

Decision:
Upload a file
Continue without uploading

SUBMIT YOUR HOMEWORK
We couldn't find that subject.
Please select the best match from the list below.

We'll send you an email right away. If it's not in your inbox, check your spam folder.

  • 1
  • 2
  • 3
Live Chats