你好，这里是arcadia13。
在这第 15 篇 MATLAB 文章中，我将解释其中一种数值计算方法，即 Runge-Kutta 方法（经典的 Runge-Kutta 方法）。
这是欧拉法的改进，与真解的误差小于欧拉法得到的近似解。
在这里，我们解释了广泛使用的“四阶龙格-库塔法”。
实际上，在 ode45() 命令中使用了 Runge-Kutta 方法。
换句话说，实现 Runge-Kutta 方法可以被认为是自己实现 ode45() 命令。
1.龙格-库塔算法
Runge-Kutta 算法与第 14 部分中解释的 Euler 方法几乎相同。所以，让我们先回顾一下欧拉的方法。
有关欧拉方法的详细信息，请参阅第 14 部分。
欧拉算法
-
确定微分方程$y_0=y(x_0)$的初始值和积分范围$a(=x_0) leq x leq b$。这为欧拉方法提供了起点 $(x,y)=(x_0, y_0)$。
-
将积分范围除以 $n$，再除以步长 $h$。即 $x_{k+1}=x_k + h
(k=0,1,...,n−1), x_0=a, x_n=b$。 -
假设$y_k=y(x_k)$，将$(x_k,y_k)$代入近似$y_{k+1}pprox y_k + hf(x_k,y_k)$得到$y_{k+1}$求
-
重复第 3 步，直到从 $y_1$ 获得 $y_n$。
Euler 方法和 Runge-Kutta 方法的区别在于近似公式。
将近似公式 $y_{k+1}pprox y_k + hf(x_k,y_k)$ 的 $f(x_k,y_k)$ 部分更改如下。
龙格-库塔法中的近似公式
egin{eqnarray}
y_{k+1} &pprox& y_k + hcfrac{k_1 + 2k_2 + 2k_3 + k_4}{6}\
k_1 &=& f(x_k, y_k)\
k_2 &=& f(x_k + rac{h}{2}, y_k + rac{h}{2}k_1)\
k_3 &=& f(x_k + rac{h}{2}, y_k + rac{h}{2}k_2)\
k_4 &=& f(x_k + h, y_k + k_3)\
end{eqnarray}
是不是处于不想看到的状态……
但是仔细看。我做的很简单，我只是通过将各种值代入$f(x,y)$来进行计算。
它与欧拉方法的作用相同，因为它被替换和计算。
实际上，欧拉法也是通过将$(x_k, y_k)$ 代入$f(x,y)$ 来计算的。我只是稍微调整了一下。
该算法的其余部分与欧拉方法完全相同。大多数 Runge-Kutta 方法都是在实现 Euler 方法时实现的。
2 Runge-Kutta 方法的实现
我们现在转向 Runge-Kutta 方法的实现。
由于我要实现 Runge-Kutta 方法（？），让我们通过结合我最喜欢的物理来解决运动方程。
主题是高中物理力学熟悉的“斜投影”。
为简单起见，我们考虑二维平面上的运动并忽略空气阻力。
让我们尝试绘制放置在平面原点的质点的轨迹，方法是在 45 度仰角的方向上发射它。
我将初始速度设置为 $10 ext{m/s}$。
% 質点の斜方投射
clear; clc;
%% 物理パラメータと運動方程式の初期条件
% 手順1,2はここで行う
x0 = 0; y0 = 0; % 原点
g = 9.8; % 重力加速度
theta = (45/180)*pi; % 45度の方向へ投げる
v0 = 10; % 初速度 m/s
% 運動方程式の初期条件
Y0 = [y0; v0*sin(theta)];
% シミュレーション時間
h = 10^(-4); %刻み幅
t = 0:h:5;
%% 運動方程式を解いてボールの軌跡を描画
x = x0 + v0*cos(theta)*t;
y = ClassicalRungeKutta(@MotionEquation, t, Y0); % 近似解
yTrue = y0 + v0*sin(theta)*t - (g*t.^2)/2; % 真の解
plot(x, y(1,:), 'LineWidth', 1.3);
hold on
plot(x, yTrue, 'g--','LineWidth', 1.3)
title('斜方投射')
legend('質点の軌跡')
xlabel('x[m]')
ylabel('y[m]')
set(gca,'FontSize',13)
ylim([0 5])
grid on
%% 斜方投射（鉛直上向き）の運動方程式
function Ydot = MotionEquation(t, Y)
g = 9.8;
nrow = length(Y);
Ydot = zeros(nrow, 1);
Ydot(1,1) = Y(2);
Ydot(2,1) = -g;
end
%% ルンゲ・クッタ法
function Yout = ClassicalRungeKutta(func, t, Y0)
% func:微分方程式の関数ハンドル
% t：積分範囲（hで分割済み）
% Y0：初期値
nrow = length(Y0);
ncol = length(t);
Y = zeros(nrow, ncol);
Y(:,1) = Y0; % 初期値の代入
h = t(2) - t(1); % 刻み幅
% 手順3,4
for k = 1:ncol-1
k1 = feval(func, t(k), Y(:,k));
k2 = feval(func, t(k) + h/2, Y(:,k) + h*k1/2);
k3 = feval(func, t(k) + h/2, Y(:,k) + h*k2/2);
k4 = feval(func, t(k) +h, Y(:,k) + h*k3);
Y(:, k+1) = Y(:, k) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
end
Yout = Y;
end
for语句的内容也可以这样写（我用@WolfMoon的代码作为参考）。
k1 = func(t(k), Y(:,k));
k2 = func(t(k) + h/2, Y(:,k) + h*k1/2);
k3 = func(t(k) + h/2, Y(:,k) + h*k2/2);
k4 = func(t(k) +h, Y(:,k) + h*k3);
Y(:, k+1) = Y(:, k) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
忽略空气阻力的斜投影有一个著名的公式（真解），可以在高中物理教科书中找到。
正好，所以比较了大概的解和官方的轨迹。
题目虽然简单，但也恰到好处。惊人。
本文主要部分到此结束。
* 鼓励深入探索龙格-库塔方法的奖励
说完Runge-Kutta法算法是这样的，我在MATLAB中实现了。
在这里，作为奖励，我将编写一个指南，指导那些有兴趣深入挖掘龙格-库塔方法的人。
用图表直观理解
无论是 Euler 方法还是 Runge-Kutta 方法，您是否曾因直觉而不确定自己在做什么？
如果您是这样的人，请访问以下网站。
这里，用图来说明近似公式的含义。视觉上吸引人的图表仍然是直观理解的最佳选择。
盯着公式看，当你真的不明白的时候，你可以通过用图表来轻松理解它。
龙格-库塔近似的推导
Runge-Kutta 方法最大的奥秘在于近似公式。那是从哪里来的？
答案是泰勒展开式。通过考虑多变量函数的泰勒展开，我们可以推导出这个近似值。
另外，虽然我一开始就写了四阶龙格-库塔法，但四阶的原因在于这个泰勒展开式，据说它是四阶龙格-库塔法，因为它甚至考虑了四阶导数项。这是。
以下是 HP 的摘要，可能会有所帮助。
概括
这一次，我们实现了 Runge-Kutta 方法，它是 Euler 方法的改进版本。
Runge-Kutta 方法非常有名，即使你自己不去实现它，知道它是知识也不会有损失。
在最后
手动求解常微分方程的方法之一是将解视为泰勒级数。
事实上，众所周知，任何单变量函数都可以用泰勒级数来表示。由于解是一个函数，我们通过认为应该将其表示为泰勒级数来求解常微分方程。
从某种意义上说，Taylor展开式可能不可避免地涉及到Runge-Kutta方法等数值计算算法。
如果您也可以阅读下一篇文章，我将非常高兴。
*如果您对本文的改进或更正有任何建议，或者您想了解类似的内容，请告诉我们。
参考网站
[1] Rai Burari，使用数值计算求解常微分方程-龙格-库塔法解释-，2020/01/19出版，2020/03/06更新，2022/10/7查看
[2] MONOist，Let's solve Differential equations using Euler and Runge-Kutta methods，2015/1/28出版，2022/10/4查看
原创声明：本文系作者授权爱码网发表，未经许可，不得转载;
原文地址：https://www.likecs.com/show-308628558.html