清空筛选条件
清空筛选条件

Runge Kutta Methods (Experimental H)

1 次查看（过去 30 天）
Mantej Sokhi
Mantej Sokhi 2022-12-3
编辑： Walter Roberson 2023-10-22
Ran in:
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)

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2022-12-4
Ran in:
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
  2 个评论
Walter Roberson
Walter Roberson 2023-10-22
Ran in:
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
Walter Roberson 2023-10-22
编辑：Walter Roberson 2023-10-22
Ran in:
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Applications 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by