Change BDF 2 step method into BDF 4 step method

%% 2step BDF method

f =@(t, y) cos(2*t + y) + (3/2)*(t - y);

y0= 1;

a = 0;

b = 1;

N0= 10;

p = 2;

j = 1;

for k = 0:5

N = N0*2^k;

NN(j) = N;

h = (b-a)/N;

hh(j) = h;

tn = a:h:b;

yn = zeros(1,length(tn));

y(1) = y0;

%%Runge-Kuta

k1=h*f(tn(1),y(1));

k2=h*f(tn(1)+h,y(1)+k1);

y(2)=y(1)+1/2*(k1+k2);

t=a;

for i=1:length(tn)-2

ynew = fzero (@(ynew) ynew - (((4/3)*y(i+1)) - ((1/3)*y(i)) + ((2/3)*f(t+2*h,ynew)*h)) , y(i));

y(i+2) = ynew;

t = t + h;

end

plot(tn, y)

hold on;

yend(j)=y(end);

j=j+1;

end

figure(1);

xlabel('t');

ylabel('y(t)');

[mk,nk]=size(yend);

for k=1:nk

end

Answer 1

[Amish]

Hi,

I see that you have a 2-Step BDF and are looking to change it to a 4-Step BDF method. In order to do so, you need to adjust the coefficients for the BDF formula and provide three additional starting values (since the 4-step BDF method requires four previous values to compute the next value). These additional starting values can be obtained using a method with sufficient order, such as the 4th-order Runge-Kutta method.

This can be done by computing the first four values of 'y' using the 4th-order Runge-Kutta method followed by a 4-step BDF method to compute the subsequent values. The coefficients (48/25, -36/25, 16/25, -3/25) and (12/25) are to be used, specific to the 4-step BDF method.

Here is a generic code snippet along the lines of your existing code:

(Make sure to adjust the rest of your code to handle the output 'yend' appropriately)

%% 4-step BDF method
f = @(t, y) cos(2*t + y) + (3/2)*(t - y);
y0 = 1;
a = 0;
b = 1;
N0 = 10;
p = 4; % 4-step BDF

j = 1;
for k = 0:5
    N = N0*2^k;
    NN(j) = N;
    h = (b-a)/N;
    hh(j) = h;
    tn = a:h:b;
    yn = zeros(1,length(tn));
    y(1) = y0;

    %% Runge-Kutta for initial values
    for i = 1:p-1
        k1 = h*f(tn(i), y(i));
        k2 = h*f(tn(i) + h/2, y(i) + k1/2);
        k3 = h*f(tn(i) + h/2, y(i) + k2/2);
        k4 = h*f(tn(i) + h, y(i) + k3);
        y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
    end

    t = a + (p-1)*h;
    for i = p:length(tn)-1
        ynew = fzero(@(ynew) ynew - (...
            (48/25)*y(i) - (36/25)*y(i-1) + (16/25)*y(i-2) - (3/25)*y(i-3) + ...
            (12/25)*h*f(t+h, ynew)), y(i)); 
        y(i+1) = ynew;
        t = t + h;
    end
    plot(tn, y)
    hold on;
    yend(j) = y(end);
    j = j+1;
end
figure(1);
xlabel('t');
ylabel('y(t)');
[mk, nk] = size(yend);
for k = 1:nk
    % Add any code here if needed for post-processing
end

Hope this helps!

[E]

Thank you! It works perfectly! I appreciate this so much! Have a nice day :)