Question

Need help with this MATLAB problem:

Using the fourth order Runge-Kutta method (KK4 to solve a first order initial value problem NOTE: This assignment is to be completed using MATLAB, and your final results including the corresponding M- iles shonma ac Given the first order initial value problem with h-time step size (i.e. ti = to + ih), then the following formula computes an approximate solution to (): i vit), where y(ti) - true value (ezact solution), (t)-f(t, v), vto) -yo (), k1 =f(t,,y.) and Forthe following segments, consider the first orderinitial value problem : 1. Solve this problem using RK4 with h 0.5, where fromt 0 to t 2 with stepsizeh-0.5, it takes 4 steps : (stepo, to 0), ti 0.5, t2-1, ts 1.5, t4 2. 2. Construct an M - file in M ATLAB to compare the above results with the ezact solution, with errors tabulated as: Numerical solution yi | Exact solution y(t.) lti Error-yi-y( 0.5 0 0.0 0.5 1.0 1.5 2.0 0.5 3. Solve the same problem using RK4 with h-0.2. 4. Solve the same problem using RK4 with h -0.05. 5. Compare the errors at t = 2 for h 05, h = 0.2 and h 0.05.

Answer 1

MATLAB Script (run it as a script, NOT from command window):

close all
clear
clc

% Parts (1) and (2)
fprintf('Parts (1) and (2)\n-----------------------------------------------------------------------------------\n')
% Exact solution
syms y(t)
eqn = diff(y,t) == y - t^2 + 1; % ODE
cond = y(0) == 0.5; % initial condition
y_exact(t) = dsolve(eqn, cond); % solving
clear y t

h = 0.5; % step size
ti = 0; tf = 2; % time span
n = (tf - ti)/h; % number of steps
t = ti:h:tf;
y(1) = 0.5; % initial condition
fprintf('%-5s\t%-30s\t%-30s\t%-20s\n', 'ti', 'Exact Solution (y(ti))', 'Numerical Solution (yi)', 'Error = |yi - y(ti)|')
for i = 1:length(t)
fprintf('%-5.2f\t%-30.8f\t%-30.8f\t%-20.16f\n', t(i), subs(y_exact, t(i)), y(i), abs(y(i) - subs(y_exact, t(i))))
k_1 = func(t(i), y(i));
k_2 = func(t(i) + 0.5*h, y(i) + 0.5*h*k_1);
k_3 = func((t(i) + 0.5*h), (y(i) + 0.5*h*k_2));
k_4 = func((t(i) + h), (y(i) + k_3*h));
y(i + 1) = y(i) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4)*h; % main update of y
end
y_end(1) = y(end-1);

% Part (3)
fprintf('\nPart (3)\n-----------------------------------------------------------------------------------\n')
h = 0.2; % step size
ti = 0; tf = 2; % time span
n = (tf - ti)/h; % number of steps
t = ti:h:tf;
y(1) = 0.5; % initial condition
fprintf('%-5s\t%-30s\t%-30s\t%-20s\n', 'ti', 'Exact Solution (y(ti))', 'Numerical Solution (yi)', 'Error = |yi - y(ti)|')
for i = 1:length(t)
fprintf('%-5.2f\t%-30.8f\t%-30.8f\t%-20.16f\n', t(i), subs(y_exact, t(i)), y(i), abs(y(i) - subs(y_exact, t(i))))
k_1 = func(t(i), y(i));
k_2 = func(t(i) + 0.5*h, y(i) + 0.5*h*k_1);
k_3 = func((t(i) + 0.5*h), (y(i) + 0.5*h*k_2));
k_4 = func((t(i) + h), (y(i) + k_3*h));
y(i + 1) = y(i) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4)*h; % main update of y
end
y_end(2) = y(end-1);

% Part (4)
fprintf('\nPart (4)\n-----------------------------------------------------------------------------------\n')
h = 0.05; % step size
ti = 0; tf = 2; % time span
n = (tf - ti)/h; % number of steps
t = ti:h:tf;
y(1) = 0.5; % initial condition
fprintf('%-5s\t%-30s\t%-30s\t%-20s\n', 'ti', 'Exact Solution (y(ti))', 'Numerical Solution (yi)', 'Error = |yi - y(ti)|')
for i = 1:length(t)
fprintf('%-5.2f\t%-30.8f\t%-30.8f\t%-20.16f\n', t(i), subs(y_exact, t(i)), y(i), abs(y(i) - subs(y_exact, t(i))))
k_1 = func(t(i), y(i));
k_2 = func(t(i) + 0.5*h, y(i) + 0.5*h*k_1);
k_3 = func((t(i) + 0.5*h), (y(i) + 0.5*h*k_2));
k_4 = func((t(i) + h), (y(i) + k_3*h));
y(i + 1) = y(i) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4)*h; % main update of y
end
y_end(3) = y(end-1);

% Part (5)
fprintf('\nPart (5)\n-----------------------------------------------------------------------------------\n')
disp('At t = 2:')
H = [0.5 0.2 0.05];
for i = 1:3
fprintf('Error (h = %.2f) = %.16f\n', H(i), abs(y_end(i) - subs(y_exact, 2)))
end

% Functions
function f = func(t,y)
f = y - t^2 + 1; % ODE
end

Output:

Parts (1) and (2)
-----------------------------------------------------------------------------------
ti    Exact Solution (y(ti))     Numerical Solution (yi)    Error = |yi - y(ti)|
0.00    0.50000000     0.50000000     0.0000000000000000
0.50    1.42563936     1.42513021     0.0005091563166026
1.00    2.64085909     2.63960266     0.0012564246376649
1.50    4.00915546     4.00681897     0.0023364947865139
2.00    5.30547195     5.30160523     0.0038667212686879

Part (3)
-----------------------------------------------------------------------------------
ti    Exact Solution (y(ti))     Numerical Solution (yi)    Error = |yi - y(ti)|
0.00    0.50000000     0.50000000     0.0000000000000000
0.20    0.82929862     0.82929333     0.0000052875865817
0.40    1.21408765     1.21407621     0.0000114405126981
0.60    1.64894060     1.64892202     0.0000185827631454
0.80    2.12722954     2.12720268     0.0000268508058225
1.00    2.64085909     2.64082269     0.0000363930417257
1.20    3.17994154     3.17989417     0.0000473683994957
1.40    3.73240002     3.73234007     0.0000599437226831
1.60    4.28348379     4.28340950     0.0000742894840376
1.80    4.81517627     4.81508569     0.0000905732140940
2.00    5.30547195     5.30536300     0.0001089498420220

Part (4)
-----------------------------------------------------------------------------------
ti    Exact Solution (y(ti))     Numerical Solution (yi)    Error = |yi - y(ti)|
0.00    0.50000000     0.50000000     0.0000000000000000
0.05    0.57686445     0.57686445     0.0000000051974046
0.10    0.65741454     0.65741453     0.0000000105939663
0.15    0.74158288     0.74158286     0.0000000161964446
0.20    0.82929862     0.82929860     0.0000000220117684
0.25    0.92048729     0.92048726     0.0000000280470360
0.30    1.01507060     1.01507056     0.0000000343095140
0.35    1.11296623     1.11296618     0.0000000408066360
0.40    1.21408765     1.21408760     0.0000000475460007
0.45    1.31834391     1.31834385     0.0000000545353690
0.50    1.42563936     1.42563930     0.0000000617826610
0.55    1.53587349     1.53587342     0.0000000692959515
0.60    1.64894060     1.64894052     0.0000000770834654
0.65    1.76472959     1.76472950     0.0000000851535717
0.70    1.88312365     1.88312355     0.0000000935147770
0.75    2.00399999     2.00399989     0.0000001021757168
0.80    2.12722954     2.12722942     0.0000001111451479
0.85    2.25267657     2.25267645     0.0000001204319358
0.90    2.38019844     2.38019831     0.0000001300450452
0.95    2.50964517     2.50964503     0.0000001399935249
1.00    2.64085909     2.64085894     0.0000001502864950
1.05    2.77367444     2.77367428     0.0000001609331278
1.10    2.90791699     2.90791682     0.0000001719426319
1.15    3.04340355     3.04340336     0.0000001833242295
1.20    3.17994154     3.17994134     0.0000001950871349
1.25    3.31732852     3.31732831     0.0000002072405288
1.30    3.45535167     3.45535145     0.0000002197935319
1.35    3.59378723     3.59378700     0.0000002327551726
1.40    3.73240002     3.73239977     0.0000002461343553
1.45    3.87094274     3.87094248     0.0000002599398232
1.50    4.00915546     4.00915519     0.0000002741801180
1.55    4.14676491     4.14676462     0.0000002888635379
1.60    4.28348379     4.28348348     0.0000003039980874
1.65    4.41901009     4.41900977     0.0000003195914271
1.70    4.55302630     4.55302597     0.0000003356508187
1.75    4.68519866     4.68519831     0.0000003521830611
1.80    4.81517627     4.81517590     0.0000003691944261
1.85    4.94259024     4.94258985     0.0000003866905849
1.90    5.06705278     5.06705237     0.0000004046765317
1.95    5.18815621     5.18815579     0.0000004231564970
2.00    5.30547195     5.30547151     0.0000004421338585

Part (5)
-----------------------------------------------------------------------------------
At t = 2:
Error (h = 0.50) = 0.0038667212686879
Error (h = 0.20) = 0.0001089498420220
Error (h = 0.05) = 0.0000004421338585