Disucssion Section Note.

Download

RKF 45 Matlab code

%% Parameters

format longE;

% define a function f
f = @(t,y) y - t^2 + 1;

% define endpoints
a = 0;
b = 2;

% tolerance
TOL = 1e-5;

% maximum and minimum step sizes
hmax = 0.25;
hmin = 0.01;

% initial condition
alpha = 0.5;

t = a;
w = alpha;
h = hmax;

%% Running the code

[t_arr, w_arr, w_hat_arr] = RKF45(f,t,w,h,b,TOL,hmax,hmin);

summary = cat(2,t_arr,w_arr,w_hat_arr,1/h*abs(w_arr-w_hat_arr));
% summary



%% RKF45 Algorithm

function [t_arr, w_arr, w_hat_arr] = RKF45(f,t,w,h,b,TOL,hmax,hmin)

    % saving the data
    t_arr = [t];
    w_arr = [w];
    w_hat_arr = [w];
    
    FLAG = 1;
    
    while FLAG == 1
        k1 = h * f(t,w);
        k2 = h * f(t+0.25*h, w+0.25*k1);
        k3 = h * f(t+3/8*h, w+3/32*k1+9/32*k2);
        k4 = h * f(t+12/13*h, w+1932/2197*k1-7200/2197*k2+7296/2197*k3);
        k5 = h * f(t+h, w+439/216*k1-8*k2+3680/513*k3-845/4104*k4);
        k6 = h * f(t+1/2*h, w-8/27*k1+2*k2-3544/2565*k3+1859/4104*k4-11/40*k5);
        
        % RK 4
        w_val     = w + 25/216*k1 + 1408/2565*k3 + 2197/4104*k4 - 1/5*k5;
        % RK 5
        w_hat_val = w + 16/135*k1 + 6656/12825*k3 + 28561/56430*k4 - 9/50*k5 + 2/55*k6;
        % Calculate R
        R = 1/h * abs(w_val - w_hat_val);
        
        if R <= TOL
            t = t+h;
            w = w_val;
            
            t_arr = [t_arr; t];
            w_arr = [w_arr; w];
            w_hat_arr = [w_hat_arr; w_hat_val];
            
            fprintf("t : %f\tw : %e\tw_hat : %e\tR : %e \n",t,w,w_hat_val,R)
        end
        
        delta = 0.84*(TOL/R)^(1/4);
        if delta <= 0.1
            h = 0.1*h;
        elseif delta >= 4
            h = 4*h;
        else
            h = delta * h;
        end
        
        if h > hmax
            h = hmax;
        end
        
        if t >= b
            FLAG = 0;
        elseif t+h > b
            h = b-t;
        elseif h < hmin
            FLAG = 0;
            fprintf("minimum h exceeded")
        end
    end
    fprintf("The number of time steps : %d\nSolution at t=%f : %f\n", length(t_arr), b, w)
end