问题描述
我正在尝试使用 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 (将#修改为@)