logo标签列表

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

matlabrunge-kutta
4

我想用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]
是的，在这个特定的例子中，f和g只是x和t的函数，但我的理解是，如果我也让它们成为y和z的函数，那么这段代码就可以应用于任何方程系统。顺便说一句，我喜欢你的编辑 :) - gelbrekt
1
“1/2M1”，“1/2M2”和“M3”是不正确的。它们不应该被包括在内，因为只有“y”在“x”的方程中被步进。 - TroyHaskin
1
问题实际上是我觉得网络上缺乏使用Rk4解决三个方程微分方程的具体示例。再次强调，我只找到了非常糟糕的具体示例和一般性描述。 - gelbrekt
1
没问题。事实上，我认为这段文字给出了一个非常糟糕的RK方法编码示例，导致代码变得极其丑陋、复杂和容易出错。如果将方程视为向量而不是单独处理每个方程，代码会更加简洁明了。但我很高兴你解决了这个问题。 - TroyHaskin
1
嗯，我不同意这会使代码变得复杂。我发现这比向量函数更直观。你能推荐一个好的页面/ PDF，用向量解决类似问题（包含具体示例）吗？ - gelbrekt
显示剩余8条评论
2个回答
6
所以你的主要问题不是正确定义 x。你使用了龙格-库塔4阶(RK4)方法传播其值，但实际上从未定义过它的导数！
在本答案底部有一个函数，可以接受任意数量的方程和它们的初始条件。这是为了满足您对三个（或更多）方程的清晰示例的需求而包含在内的。
供参考，这些方程可以直接从标准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函数给出了相同的结果。

results comparison

函数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
1

好的，事实证明只是一个小错误，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]
1
我想指出，为了使这个工作正常运行，不需要定义e；使用解向量的原始值也是可行的。然而，定义e会显著降低错误的可能性，因为实际的函数形式被隐藏在函数调用之下。 - TroyHaskin
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接