实现龙格库塔系统的DE

Question

我试图在MATLAB中实现Runge-Kutta Method for Systems of DEs。我没有得到correct answers，我不确定代码中是否有错误或者我用来运行它的命令。

这里是我的代码：

function RKSystems(a, b, m, N, alpha, f) 
    h = (b - a)/N; 
    t = a; 
    w = zeros(1, m); 

    for j = 1:m 
     w(j) = alpha(j); 
    end 

    fprintf('t = %.2f;', t); 
    for i = 1:m 
     fprintf(' w(%d) = %.10f;', i, w(i)); 
    end 
    fprintf('\n'); 

    k = zeros(4, m); 
    for i = 1:N 
     for j = 1:m 
      k(1, j) = h*f{j}(t, w); 
     end 

     for j = 1:m 
      k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1)); 
     end 

     for j = 1:m 
      k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2)); 
     end 

     for j = 1:m 
      k(4, j) = h*f{j}(t + h, w + k(3)); 
     end 

     for j = 1:m 
      w(j) = w(j) + (k(1, j) + 2*k(2, j) + 2*k(3, j) + k(4, j))/6; 
     end 

     t = a + i*h; 

     fprintf('t = %.2f;', t); 
     for k = 1:m 
      fprintf(' w(%d) = %.10f;', k, w(k)); 
     end 
     fprintf('\n'); 

    end 
end 

我想测试一下在this problem。 这是我的命令和输出：

>> U1 = @（T，U）3 * U（1）+ 2 * U（2） - （2 * T^2 + 1）* EXP（ 2 * T）; （1）+ u（2）+（t^2 + 2 * t-4）* exp（2 * t）;其中，U2 =

>> a = 0; b = 1; alpha = [1 1]; m = 2; h = 0.2; N =（b-a）/ h; RKSystems（a，b，m，N，alpha，{U1 U2}）;

t = 0.00; w（1）= 1.0000000000; w（2）= 1.0000000000;

t = 0.20; w（1）= 2.2930309680; w（2）= 1.6186020410;

t = 0.40; w（1）= 5.0379629523; w（2）= 3.7300162165;

t = 0.60; w（1）= 11.4076339762; w（2）= 9.7009491301;

t = 0.80; w（1）= 27.0898576892; w（2）= 25.6603894354;

t = 1.00; w（1）= 67.1832886708; w（2）= 67.6103125539;

Answer 1

我的问题是在下面的代码行：

k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1));

k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2));

k(4, j) = h*f{j}(t + h, w + k(3));

我期待k(1)由米添加的k（4整个第一行矩阵）到矩阵w（一个1乘m矩阵），但它只是添加第一行的第一个元素。为了解决这个问题，我修改了行，如下所示：

k(2, j) = h*f{j}(t + h/2, w + (1/2)*k(1, :));

k(3, j) = h*f{j}(t + h/2, w + (1/2)*k(2, :));

k(4, j) = h*f{j}(t + h, w + k(3, :));

Answer 2

所以...这是我该怎么做，我很难阅读代码片段中发生了什么，但我希望这可以帮助你。此外，matlab内置了rk解算器，我建议熟悉这些解算器。

%rk4 solver 
dt = .2; 
t = 0:dt:1; 
u = zeros(2,numel(t)); 
u(:,1) = 1; 

for jj = 2:numel(t) 
    u_ = u(:,jj-1); 
    t_ = t(jj-1); 
    fa = rhs(u_,t_); 
    fb = rhs(u_+dt/2.*fa,t_+dt/2); 
    fc = rhs(u_+dt/2.*fb,t_+dt/2); 
    fd = rhs(u_+dt.*fc,t_+dt); 
    u(:,jj) = u(:,jj-1) + dt/6*(fa+2*fb+2*fc+fd); 
end 
disp([t;u]') 

和rhs.m如下：

function dudt = rhs(u,t) 
dudt = [(3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t)); 
     (4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t))]; 

这将返回以下：或者

>> rkorder4 
    0 1.0000 1.0000 
0.2000 2.1204 1.5070 
0.4000 4.4412 3.2422 
0.6000 9.7391 8.1634 
0.8000 22.6766 21.3435 
1.0000 55.6612 56.0305 

，你可以调用ODE45的东西，如：

dt = .2; 
t = 0:dt:1; 
rhs=@(t,u)[(3*u(1) + 2*u(2) - (2*t^2 + 1)*exp(2*t)); 
     (4*u(1) + u(2) + (t^2 + 2*t - 4)*exp(2*t))]; 

[ts,us]=ode45(@(t,u) rhs(t,u),t,[1 1],[]); 
disp([ts us]); 

它给你：

    0 1.000000000000000 1.000000000000000 
    0.200000000000000 2.125018862140859 1.511597928697035 
    0.400000000000000 4.465156492097592 3.266022178547346 
    0.600000000000000 9.832481410895053 8.256418221678395 
    0.800000000000000 23.003033457636558 21.669270713784108 
    1.000000000000000 56.738351357036301 57.106230777717208 

这是非常接近你从我的代码得到的。两者之间的协议可以通过减少时间步长dt来增加。在t的低值时，它们将始终保持最大一致性，两者之间的差异随着t的增加而增加。这两种实现都与分析解决方案非常接近。