Runge Kutta Methods (Experimental H)

I am trying to see the differences in my plot with different values of h. I would like to just run the program once with different values of h and have all the graphs in the same plot . Do I just put all my values of h in an array ? Below is my code:

h = 0.005;
tfinal = 55;

y(1) = 1;
t(1) = 1;

f = @(t,y) -y^2;

for i = 1:ceil(tfinal/h)

    t(i+1) = t(i) + h;
    k1 = f(t(i),y(i));
    k2 = f(t(i)+0.5*h,y(i)+0.5*k1*h);
    k3 = f(t(i)+0.5*h,y(i)+0.5*k2*h);
    y(i+1) = y(i) + h/6* (k1 + 2*k2 + 2*k3);
    
end

plot(t,y)
xlabel('t')
ylabel('y')
set(gca,'Fontsize',16)

Answer 1

[Walter Roberson]

Indistinguishable results on the scale of h that I experimented with

h = 0.002:0.001:0.007;
num_h = length(h);

tfinal = 55;

y = ones(1, num_h);
t = ones(1, num_h);

f = @(t,y) -y.^2;

for i = 1 : max(ceil(tfinal./h))

    t(i+1, :) = t(i,:) + h;
    k1 = f(t(i,:),y(i,:));
    k2 = f(t(i,:)+0.5*h, y(i,:)+0.5*k1.*h);
    k3 = f(t(i,:)+0.5*h, y(i,:)+0.5*k2.*h);
    y(i+1,:) = y(i,:) + h/6 .* (k1 + 2*k2 + 2*k3);
end

for K = 1 : num_h
    subplot(num_h,1,K)
    plot(t(:,K), y(:,K))
    title(string(h(K)));
    xlabel('t')
    ylabel('y')
    set(gca,'Fontsize',16)
end

[Walter Roberson]

You asked for them all on one plot, so here goes.

It looks like only a single plot because the values are so close together.

h = 0.002:0.001:0.007;
num_h = length(h);

tfinal = 55;

y = ones(1, num_h);
t = ones(1, num_h);

f = @(t,y) -y.^2;

for i = 1 : max(ceil(tfinal./h))

    t(i+1, :) = t(i,:) + h;
    k1 = f(t(i,:),y(i,:));
    k2 = f(t(i,:)+0.5*h, y(i,:)+0.5*k1.*h);
    k3 = f(t(i,:)+0.5*h, y(i,:)+0.5*k2.*h);
    y(i+1,:) = y(i,:) + h/6 .* (k1 + 2*k2 + 2*k3);
end

for K = 1 : num_h
    plot(t(:,K), y(:,K), 'DisplayName', string(h(K)));
    hold on
end
xlabel('t')
ylabel('y')
set(gca,'Fontsize',16)
legend show

[Walter Roberson]

The plots are not all the same -- they just are so close that when plotted together you cannot distinguish them. Below what is plotted is the difference between the output and the first output -- but because they are all operating on different time vectors first you have to interpolate them to the same time vector for comparison.

The difference is at most 3e-5 which is not something that is going to show up on a plot with a y range of 0 to 1 like your original plot.

h = 0.002:0.001:0.007;
num_h = length(h);

tfinal = 55;

y = ones(1, num_h);
t = ones(1, num_h);

f = @(t,y) -y.^2;

for i = 1 : max(ceil(tfinal./h))

    t(i+1, :) = t(i,:) + h;
    k1 = f(t(i,:),y(i,:));
    k2 = f(t(i,:)+0.5*h, y(i,:)+0.5*k1.*h);
    k3 = f(t(i,:)+0.5*h, y(i,:)+0.5*k2.*h);
    y(i+1,:) = y(i,:) + h/6 .* (k1 + 2*k2 + 2*k3);
end

for K = 2 : num_h
    Y = interp1(t(:,K), y(:,K), t(:,1));
    plot(t(:,1), Y - y(:,1), 'DisplayName', string(h(K)));
    hold on
end
xlabel('t')
ylabel('y')
set(gca,'Fontsize',16)
legend show