使用 Runge Kutta 4 求解方程组：Matlab

Question

我想在 Matlab 中使用 Runge Kutta 4 方法求解三个微分方程组（不允许使用 Ode45）。

经过长时间的查找，我在网上找到的要么是难以理解的示例，要么是根本不包含示例的一般解释。我想要一个关于如何正确实施我的解决方案的具体示例，或者我可以构建的类似问题的解决方案。

我已经走了很远； 我当前的代码在大多数组件上输出了一个包含 2 个正确小数的矩阵，我对此非常满意。

但是，当步长减小时，误差会变得很大。我知道我创建的 for 循环并不完全正确。我可能错误地定义了函数，但我很确定如果对 for 循环进行一些小的更改，问题就解决了，因为它似乎已经很好地解决了当前状态下的方程组。

clear all, close all, clc

%{
____________________TASK:______________________
Solve the system of differential equations below 
in the interval 0<t<1, with stepsize h = 0.1.
x'= y                 x(0)=1
y'= -x-2e^t+1         y(0)=0   ,   where x=x(t), y=y(t),  z=z(t)  
z'= -x - e^t + 1      z(0)=1

THE EXACT SOLUTIONS for x y and z can be found in this pdf:
archives.math.utk.edu/ICTCM/VOL16/C029/paper.pdf
_______________________________________________

%}

h = 0.1;
t    = 0:h:1
N = length(t);

%Defining the functions
x    = zeros(N,1);%I am not entierly sure if x y z are supposed to be defined in this way.
y    = zeros(N,1)
z    = zeros(N,1)

f = @(t, x, y, z) -x-2*exp(t)+1;%Question: Do i need a function for x here as well??
g = @(t, x, y, z) -x - exp(t) + 1;

%Starting conditions
x(1) = 1; 
y(1) = 0;
z(1) = 1;

for i = 1:(N-1)
    K1     = h * ( y(i));%____I think z(i) is supposed to be here, but i dont know in what way. 
    L1     = h * f( t(i)          , x(i)        , y(i) ,      z(i));
    M1     = h * g( t(i)          , x(i)        , y(i) ,      z(i));

    K2     = h *  (y(i) + 1/2*L1 + 1/2*M1);%____Again, z(i) should probably be here somewhere. 
    L2     = h * f(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    M2     = h * g(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);

    K3     = h *  (y(i) + 1/2*L2 + 1/2*M2);%____z(i). Should it just be added, like "+z(i)" ? 
    L3     = h * f(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    M3     = h * g(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);

    K4     = h *  (y(i) + L3 + M3);%_____z(i) ... ? 
    L4     = h * f( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    M4     = h * g( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);

    x(i+1) = x(i)+1/6*(K1+2*K2+2*K3+K4);                                                                             
    y(i+1) = y(i)+1/6*(L1+2*L2+2*L3+L4);
    z(i+1) = z(i)+1/6*(M1+2*M2+2*M3+M4);
end

Answer_Matrix = [t' x y z]

Answer 1

所以你的主要问题是没有正确定义x。您使用 Runge Kutta 4 (RK4) 方法传播它的价值，但从未真正定义它的导数是什么！

---

此答案的底部是一个函数，它可以采用任何给定数量的方程及其初始条件。 包含此函数是为了满足您对三个（或更多）的清晰示例的需求方程。

作为参考，可以直接从here描述的标准RK4方法中提取方程。

---

工作脚本

这与您的类似，但使用稍微清晰的命名约定和结构。

% Initialise step-size variables
h = 0.1;
t = (0:h:1)';
N = length(t);

% Initialise vectors
x = zeros(N,1);    y = zeros(N,1);    z = zeros(N,1);
% Starting conditions
x(1) = 1;     y(1) = 0;    z(1) = 1;

% Initialise derivative functions
dx = @(t, x, y, z) y;                  % dx = x' = dx/dt
dy = @(t, x, y, z) - x -2*exp(t) + 1;  % dy = y' = dy/dt
dz = @(t, x, y, z) - x -  exp(t) + 1;  % dz = z' = dz/dt

% Initialise K vectors
kx = zeros(1,4); % to store K values for x
ky = zeros(1,4); % to store K values for y
kz = zeros(1,4); % to store K values for z
b = [1 2 2 1];   % RK4 coefficients

% Iterate, computing each K value in turn, then the i+1 step values
for i = 1:(N-1)        
    kx(1) = dx(t(i), x(i), y(i), z(i));
    ky(1) = dy(t(i), x(i), y(i), z(i));
    kz(1) = dz(t(i), x(i), y(i), z(i));

    kx(2) = dx(t(i) + (h/2), x(i) + (h/2)*kx(1), y(i) + (h/2)*ky(1), z(i) + (h/2)*kz(1));
    ky(2) = dy(t(i) + (h/2), x(i) + (h/2)*kx(1), y(i) + (h/2)*ky(1), z(i) + (h/2)*kz(1));
    kz(2) = dz(t(i) + (h/2), x(i) + (h/2)*kx(1), y(i) + (h/2)*ky(1), z(i) + (h/2)*kz(1));

    kx(3) = dx(t(i) + (h/2), x(i) + (h/2)*kx(2), y(i) + (h/2)*ky(2), z(i) + (h/2)*kz(2));
    ky(3) = dy(t(i) + (h/2), x(i) + (h/2)*kx(2), y(i) + (h/2)*ky(2), z(i) + (h/2)*kz(2));
    kz(3) = dz(t(i) + (h/2), x(i) + (h/2)*kx(2), y(i) + (h/2)*ky(2), z(i) + (h/2)*kz(2));

    kx(4) = dx(t(i) + h, x(i) + h*kx(3), y(i) + h*ky(3), z(i) + h*kz(3));
    ky(4) = dy(t(i) + h, x(i) + h*kx(3), y(i) + h*ky(3), z(i) + h*kz(3));
    kz(4) = dz(t(i) + h, x(i) + h*kx(3), y(i) + h*ky(3), z(i) + h*kz(3));

    x(i+1) = x(i) + (h/6)*sum(b.*kx);       
    y(i+1) = y(i) + (h/6)*sum(b.*ky);       
    z(i+1) = z(i) + (h/6)*sum(b.*kz);        
end    

% Group together in one solution matrix
txyz = [t,x,y,z];

---

作为函数实现

您想要可以“应用于任何方程系统”的代码。为了使您的脚本更有用，让我们利用向量输入，其中每个变量都在自己的行上，然后将其变成一个函数。结果与 Matlab 自己的ode45 相当（在其调用方式上）。

% setup
odefun = @(t, y) [y(2); -y(1) - 2*exp(t) + 1; -y(1) - exp(t) + 1];
y0 = [1;0;1];
% ODE45 solution
[T, Y] = ode45(odefun, [0,1], y0);
% Custom RK4 solution
t = 0:0.1:1;
y = RK4(odefun, t, y0);
% Compare results
figure; hold on;
plot(T, Y); plot(t, y, '--', 'linewidth', 2)

您可以看到RK4 函数（如下）给出了与ode45 函数相同的结果。

函数 RK4 只是上述脚本的“浓缩”版本，它适用于您想使用的任何方程。为了广泛使用，您可能希望在函数中包含输入检查。为了清楚起见，我忽略了这一点。

function y = RK4(odefun, tspan, y0)
% ODEFUN contains the ode functions of the system
% TSPAN  is a 1D vector of equally spaced t values
% Y0     contains the intial conditions for the system variables

    % Initialise step-size variables
    t = tspan(:); % ensure column vector = (0:h:1)';
    h = t(2)-t(1);% define h from t
    N = length(t);

    % Initialise y vector, with a column for each equation in odefun
    y = zeros(N, numel(y0));
    % Starting conditions
    y(1, :) = y0(:)';  % Set intial conditions using row vector of y0

    k = zeros(4, numel(y0));              % Initialise K vectors
    b = repmat([1 2 2 1]', 1, numel(y0)); % RK4 coefficients

    % Iterate, computing each K value in turn, then the i+1 step values
    for i = 1:(N-1)        
        k(1, :) = odefun(t(i), y(i,:));        
        k(2, :) = odefun(t(i) + (h/2), y(i,:) + (h/2)*k(1,:));        
        k(3, :) = odefun(t(i) + (h/2), y(i,:) + (h/2)*k(2,:));        
        k(4, :) = odefun(t(i) + h, y(i,:) + h*k(3,:));

        y(i+1, :) = y(i, :) + (h/6)*sum(b.*k);    
    end    
end

Answer 2

好吧，事实证明这只是一个小错误，其中 x 变量未定义为 y 的函数（根据问题，x'(t)=y。

所以：下面是一个如何在 matlab 中使用 Runge Kutta 4 求解微分方程组的具体示例：

clear all, close all, clc

%{
____________________TASK:______________________
Solve the system of differential equations below 
in the interval 0<t<1, with stepsize h = 0.1.
x'= y                 x(0)=1
y'= -x-2e^t+1         y(0)=0   ,   where x=x(t), y=y(t),  z=z(t)  
z'= -x - e^t + 1      z(0)=1

THE EXACT SOLUTIONS for x y and z can be found in this pdf:
archives.math.utk.edu/ICTCM/VOL16/C029/paper.pdf
_______________________________________________

%}

%Step-size
h = 0.1;
t    = 0:h:1
N = length(t);

%Defining the vectors where the answer is stored. 
x    = zeros(N,1);
y    = zeros(N,1)
z    = zeros(N,1)

%Defining the functions
e = @(t, x, y, z) y;
f = @(t, x, y, z) -x-2*exp(t)+1;
g = @(t, x, y, z) -x - exp(t) + 1;

%Starting/initial conditions
x(1) = 1; 
y(1) = 0;
z(1) = 1;

for i = 1:(N-1)
    K1     = h * e( t(i)          , x(i)        , y(i) ,      z(i));
    L1     = h * f( t(i)          , x(i)        , y(i) ,      z(i));
    M1     = h * g( t(i)          , x(i)        , y(i) ,      z(i));

    K2     = h * e(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    L2     = h * f(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);
    M2     = h * g(t(i) + 1/2*h,  x(i)+1/2*K1 , y(i)+1/2*L1 , z(i)+1/2*M1);

    K3     = h * e(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    L3     = h * f(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);
    M3     = h * g(t(i) + 1/2*h,  x(i) + 1/2*K2 , y(i) + 1/2*L2 , z(i) + 1/2*M2);

    K4     = h * e( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    L4     = h * f( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);
    M4     = h * g( t(i)+h    ,  x(i)+K3     , y(i)+L3,     z(i)+M3);

    x(i+1) = x(i)+1/6*(K1+2*K2+2*K3+K4);                                                                             
    y(i+1) = y(i)+1/6*(L1+2*L2+2*L3+L4);
    z(i+1) = z(i)+1/6*(M1+2*M2+2*M3+M4);
end

Answer_Matrix = [t' x y z]