Ode-system 输出在 Matlab 中给出空值

编程问答 2022-05-31

问题描述

我正在尝试使用 Runga Katta（四阶）方法和我选择的步骤计算这个微分方程组的数值解（在区间 [0,2] 上）。

% This is not just equations,not MATLAB code
dy1/dt = 1.3 * (y2-y1) + 1.04 * 10^4 * k * y2
dy2/dt = 1.88 * 10^3 * (y4-(1+k))
dy3/dt = 1752+(266.1 * y2)-(269.3 * y3)
dy4/dt = 0.1+(320 * y2)-(321 * y4)
k = 6 * 10^(-4) * exp(20.7 - 15 * 10^3/y1)
y0 = transpose([759.167,600,0.1])

我采用了来自 this answer 的以下代码，该代码运行流畅但输出零表（除非单独输入）。我尝试了不同的步长，但无济于事。

clc

h = 0.2; % Step size.
t = (0:h:2)'; 
N = length(t);

y1 = zeros(N,1);    
y2 = zeros(N,1);    
y3 = zeros(N,1);
y4 = zeros(N,1);

% Initial conditions
y1(1) = 759.167;     
y2(1) = 0;    
y3(1) = 600;
y4(1) = 0.1;

% Initialise derivative functions
dy1 = @(t,y1,y2,y3,y4) (1.3*(y2-y1))+1.04*10^4*6*10^(-4)*exp(20.7 - (15*10^(3)/y1))*y2;
dy2 = @(t,y4) 1.88*10^3*(y4-(1 + 6*10^(-4)*exp(20.7 - (15*10^(3)/y1))));
dy3 = @(t,y4) 1752+(266.1*y2)-(269.3*y3);
dy4 = @(t,y4) 0.1+(320*y2)-(321*y4);

% Initialise K vectors
ky1 = zeros(1,4);
ky2 = zeros(1,4);
ky3 = zeros(1,4);
ky4 = zeros(1,4);

% RK4 coefficients
b = [1 2 2 1];

% Iterate,computing each K value in turn,then the i+1 step values
for i = 1:(N-1)        
    ky1(1) = dy1(t(i),y1(i),y2(i),y3(i),y4(i));
    ky2(1) = dy2(t(i),y4(i));
    ky3(1) = dy3(t(i),y4(i));
    ky4(1) = dy4(t(i),y4(i));

    ky1(2) = dy1(t(i) + (h/2),y1(i) + (h/2)*ky1(1),y2(i) + (h/2)*ky2(1),y3(i) + (h/2)*ky3(1),y4(i) + (h/2)*ky4(1));
    ky2(2) = dy2(t(i) + (h/2),y4(i) + (h/2)*ky4(1));
    ky3(2) = dy3(t(i) + (h/2),y4(i) + (h/2)*ky4(1));
    ky4(2) = dy4(t(i) + (h/2),y4(i) + (h/2)*ky4(1));

    ky1(3) = dy1(t(i) + (h/2),y1(i) + (h/2)*ky1(2),y2(i) + (h/2)*ky2(2),y3(i) + (h/2)*ky3(2),y4(i) + (h/2)*ky4(2));
    ky2(3) = dy2(t(i) + (h/2),y4(i) + (h/2)*ky4(2));
    ky3(3) = dy3(t(i) + (h/2),y4(i) + (h/2)*ky4(2));
    ky4(3) = dy4(t(i) + (h/2),y4(i) + (h/2)*ky4(2));

    ky1(4) = dy1(t(i) + h,y1(i) + h*ky1(3),y2(i) + h*ky2(3),y3(i) + h*ky3(3),y4(i) + h*ky4(3));
    ky2(4) = dy2(t(i) + h,y4(i) + h*ky4(3));
    ky3(4) = dy3(t(i) + h,y4(i) + h*ky4(3));
    ky4(4) = dy4(t(i) + h,y4(i) + h*ky4(3));

    %y1(i+1)
    y1(i+1) = y1(i) + (h/6)*sum(b.*ky1);       
    y2(i+1) = y2(i) + (h/6)*sum(b.*ky2);       
    y3(i+1) = y3(i) + (h/6)*sum(b.*ky3);
    y4(i+1) = y4(i) + (h/6)*sum(b.*ky4);
end

txyz = [t,y3]

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

matlabmatlabnumerical-methodsrunge-kutta