使用 Runge Kutta 4 求解方程组:Matlab [英] Solve a system of equations with Runge Kutta 4: Matlab
问题描述
我想在 Matlab 中使用 Runge Kutta 4 方法求解三个微分方程组(不允许使用 Ode45).
经过长时间的查找，我在网上找到的要么是无法理解的示例，要么是根本不包含示例的一般解释.我想要一个关于如何正确实施我的解决方案的具体例子，或者我可以建立的类似问题的解决方案.
我已经走了很远；我当前的代码在大多数组件上输出了一个包含 2 个正确小数的矩阵，我对此非常满意.
然而，当步长减小时，误差会变得很大.我知道我创建的 for 循环并不完全正确.我可能错误地定义了函数，但我很确定如果对 for 循环进行一些细微的更改，问题就会得到解决，因为它似乎已经在当前状态下很好地解决了方程系统.
全部清除，全部关闭，clc%{____________________任务:______________________求解下面的微分方程组在0＜t＜1区间内，步长h＝0.1.x'= y x(0)=1y'= -x-2e^t+1 y(0)=0 , 其中 x=x(t), y=y(t), z=z(t)z'= -x - e^t + 1 z(0)=1可以在此 pdf 中找到 x y 和 z 的确切解决方案:archives.math.utk.edu/ICTCM/VOL16/C029/paper.pdf_________________________________________________________%}h = 0.1；t = 0:h:1N = 长度(t)；%定义函数x = zeros(N,1);%我不确定 x y z 是否应该以这种方式定义.y = 零点(N,1)z = 零(N,1)f = @(t, x, y, z) -x-2*exp(t)+1;%问题:我是否也需要 x 的函数?g = @(t, x, y, z) -x - exp(t) + 1;%起始条件x(1) = 1;y(1) = 0;z(1) = 1;对于 i = 1:(N-1)K1 = h * ( 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 * (y(i) + 1/2*L1 + 1/2*M1);%____同样，z(i) 应该在某个地方.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).是否应该添加，例如 "+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)；结尾Answer_Matrix = [t' x y z]所以你的主要问题是没有正确定义 x.您使用 Runge Kutta 4 (RK4) 方法传播其值，但从未真正定义其导数是什么！
在这个答案的底部是一个函数，它可以采用任何给定数量的方程及其初始条件. 包含这个函数是为了满足您对三个(或更多)的清晰示例的需求等式.
作为参考，方程可以直接从描述的标准 RK4 方法中提取
RK4 函数只是上述脚本的精简"版本，它适用于您想要使用的许多方程.对于广泛使用，您可能希望在函数中包含输入检查.为清楚起见，我已将其省略.
function y = RK4(odefun, tspan, y0)% ODEFUN 包含系统的ode函数% TSPAN 是等距 t 值的一维向量% Y0 包含系统变量的初始条件% 初始化步长变量t = tspan(:);% 确保列向量 = (0:h:1)';h = t(2)-t(1);% 从 t 定义 hN = 长度(t)；% 初始化 y 向量，odefun 中的每个方程都有一列y = zeros(N, numel(y0));% 起始条件y(1, :) = y0(:)';% 使用 y0 的行向量设置初始条件k = zeros(4, numel(y0));% 初始化 K 个向量b = repmat([1 2 2 1]', 1, numel(y0));% RK4 系数% 迭代，依次计算每个K值，然后i+1步值对于 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);结尾结尾I want to solve a system of THREE differential equations with the Runge Kutta 4 method in Matlab (Ode45 is not permitted).
After a long time spent looking, all I have been able to find online are either unintelligible examples or general explanations that do not include examples at all. I would like a concrete example on how to implement my solution properly, or the solution to a comparable problem which I can build on.
I have come quite far; my current code spits out a matrix with 2 correct decimals on most of the components, which I am quite happy with.
However, when the step-size is decreased, the errors become enormous. I know the for-loop I have created is not entirely correct. I may have defined the functions incorrectly, but I am quite certain that the problem is solved if some minor changes are made to the for-loop because it appears to be solving the equation-system fairly well already in its current state.
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]
So your main issue was not defining x properly. You were propagating its value using the Runge Kutta 4 (RK4) method, but never actually defined what its derivative was!
At the bottom of this answer is a function which can take any given number of equations and their initial conditions. This has been included to address your need for a clear example for three (or more) equations.
For reference, the equations can be directly lifted from the standard RK4 method described here.
Working Script
This is comparable to yours, but uses slightly clearer naming conventions and structure.
% 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];
Implemented as function
You wanted code which can "be applied to any equation system". To make your script more usable, let's take advantage of vector inputs, where each variable is on its own row, and then make it into a function. The result is something comparable (in how it is called) to Matlab's own 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)
You can see that the RK4 function (below) gives the same result of the ode45 function.
The function RK4 is simply a "condensed" version of the above script, it will work for however many equations you want to use. For broad use, you would want to include input-checking in the function. I have left this out for clarity.
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
这篇关于使用 Runge Kutta 4 求解方程组:Matlab的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
