Блог

Я хочу создать несколько графиков для решения ODE, чтобы изучить влияние размера шага на сходимость, как я могу настроить следующий код MATLAB, чтобы иметь несколько графиков для нескольких значений ‘h’ на одном графике?

 x(1)=0;
y(1)=2;
z(1)=4;
h=0.2;

for i=1:5
    x(i 1)=i*h;
    ky1=-2*y(i) 5*exp(-x(i));
    kz1=-(1/3)*y(i)*z(i)^2;
    ky2=-2*(y(i) ky1*h/2) 5*(exp(-(x(i) h/2)));
    kz2=-(1/3)*(y(i) ky1*h/2)*(z(i) kz1*h/2)^2;
    ky3=-2*(y(i) ky2*h/2) 5*(exp(-(x(i) h/2)));
    kz3=-(1/3)*(y(i) ky2*h/2)*(z(i) kz2*h/2)^2;
    ky4=-2*(y(i) ky3*h) 5*(exp(-(x(i) h)));
    kz4=-(1/3)*(y(i) ky3*h)*(z(i) kz3*h)^2;
    y(i 1)=y(i) (ky1 2*ky2 2*ky3 ky4)*h/6;
    z(i 1)=z(i) (kz1 2*kz2 2*kz3 kz4)*h/6;
end

step = [0:5]';
x_i=x';
y_i=y';
z_i=z';

table(step, x_i, y_i, z_i)

plot(x,y,'--mo');
hold on
plot(x,z,'--ro');

legend('y(x) 4th-order RK method (h=0.2)','z(x) 4th-order RK method (h=0.2)');
xlabel('x')
ylabel('y amp; z')
grid on

 

Вы можете добавить внешний цикл for, который проходит через каждый индекс h = [1,0.5,0.2,0.1,0.01] . Здесь я помещаю весь код, который выполняет вычисления, в вызываемую функцию Calculate() . Циклическое выполнение этой функции с использованием значений в векторе / матрице h позволит вам получить графики относительно h .

Полный сценарий: все на одном графике

 x(1)=0;
y(1)=2;
z(1)=4;
h = [1,0.5,0.2,0.1,0.01];

for Index = 1: length(h)
[x,y,z] = Calculate(x,y,z,h(Index));

plot(x,y,'--mo');
hold on
plot(x,z,'--ro');

end

legend('y(x) 4th-order RK method (h=0.2)','z(x) 4th-order RK method (h=0.2)');
xlabel('x');
ylabel('y amp; z');
grid on


function [x,y,z] = Calculate(x,y,z,h)

for i=1:5
    x(i 1)=i*h;
    ky1=-2*y(i) 5*exp(-x(i));
    kz1=-(1/3)*y(i)*z(i)^2;
    ky2=-2*(y(i) ky1*h/2) 5*(exp(-(x(i) h/2)));
    kz2=-(1/3)*(y(i) ky1*h/2)*(z(i) kz1*h/2)^2;
    ky3=-2*(y(i) ky2*h/2) 5*(exp(-(x(i) h/2)));
    kz3=-(1/3)*(y(i) ky2*h/2)*(z(i) kz2*h/2)^2;
    ky4=-2*(y(i) ky3*h) 5*(exp(-(x(i) h)));
    kz4=-(1/3)*(y(i) ky3*h)*(z(i) kz3*h)^2;
    y(i 1)=y(i) (ky1 2*ky2 2*ky3 ky4)*h/6;
    z(i 1)=z(i) (kz1 2*kz2 2*kz3 kz4)*h/6;
end

step = [0:5]';
x_i=x';
y_i=y';
z_i=z';

table(step, x_i, y_i, z_i)

end
 

Полный сценарий: в подзаголовках

 x(1)=0;
y(1)=2;
z(1)=4;
h = [1,0.5,0.2,0.1,0.01];

for Index = 1: length(h)
[x,y,z] = Calculate(x,y,z,h(Index));

subplot(length(h),1,Index); plot(x,y,'--mo');
hold on
subplot(length(h),1,Index); plot(x,z,'--ro');
legend(['y(x) 4th-order RK method (h=' num2str(h(Index)) ')'],['z(x) 4th-order RK method (h=' num2str(h(Index)) ')']);
xlabel('x');
ylabel('y amp; z');
grid on

end

function [x,y,z] = Calculate(x,y,z,h)

for i=1:5
    x(i 1)=i*h;
    ky1=-2*y(i) 5*exp(-x(i));
    kz1=-(1/3)*y(i)*z(i)^2;
    ky2=-2*(y(i) ky1*h/2) 5*(exp(-(x(i) h/2)));
    kz2=-(1/3)*(y(i) ky1*h/2)*(z(i) kz1*h/2)^2;
    ky3=-2*(y(i) ky2*h/2) 5*(exp(-(x(i) h/2)));
    kz3=-(1/3)*(y(i) ky2*h/2)*(z(i) kz2*h/2)^2;
    ky4=-2*(y(i) ky3*h) 5*(exp(-(x(i) h)));
    kz4=-(1/3)*(y(i) ky3*h)*(z(i) kz3*h)^2;
    y(i 1)=y(i) (ky1 2*ky2 2*ky3 ky4)*h/6;
    z(i 1)=z(i) (kz1 2*kz2 2*kz3 kz4)*h/6;
end

step = [0:5]';
x_i=x';
y_i=y';
z_i=z';

table(step, x_i, y_i, z_i)

end
 

Запускался с использованием MATLAB R2019b