用 Runge Kutta 4 求解方程组：Matlab
Solve a system of equations with Runge Kutta 4: Matlab
我想用 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]
好的，事实证明这只是一个小错误，其中 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]
所以您的主要问题是没有正确定义 x。您使用 Runge Kutta 4 (RK4) 方法传播它的值，但从未真正定义它的导数是什么！
在这个答案的底部是一个函数，它可以采用任意给定数量的方程及其初始条件。包含这个是为了满足您对一个清晰示例的需求三个（或更多）方程式。
作为参考，可以直接从描述的标准 RK4 方法中提取方程式 here。
工作脚本
这与您的相当，但使用了更清晰的命名约定和结构。
% 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];
作为函数实现
您想要的代码可以 "be applied to any equation system"。为了使您的脚本更有用，让我们利用向量输入，其中每个变量都在其自己的行中，然后将其变成一个函数。结果与 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 只是上述脚本的 "condensed" 版本，它适用于任何您想使用的方程式。为了广泛使用，您可能希望在函数中包含输入检查。为了清楚起见，我将其省略。
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
我想用 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]
好的，事实证明这只是一个小错误，其中 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]
所以您的主要问题是没有正确定义 x。您使用 Runge Kutta 4 (RK4) 方法传播它的值，但从未真正定义它的导数是什么！
在这个答案的底部是一个函数，它可以采用任意给定数量的方程及其初始条件。包含这个是为了满足您对一个清晰示例的需求三个（或更多）方程式。
作为参考，可以直接从描述的标准 RK4 方法中提取方程式 here。
工作脚本
这与您的相当，但使用了更清晰的命名约定和结构。
% 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];
作为函数实现
您想要的代码可以 "be applied to any equation system"。为了使您的脚本更有用，让我们利用向量输入，其中每个变量都在其自己的行中，然后将其变成一个函数。结果与 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 只是上述脚本的 "condensed" 版本，它适用于任何您想使用的方程式。为了广泛使用，您可能希望在函数中包含输入检查。为了清楚起见，我将其省略。
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