龙格库塔法的原理

在百度百科中是这么解释的：在各种龙格－库塔法当中有一个方法十分常用，以至于经常被称为“RK4”或者就是“龙格－库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。

令初值问题表述如下。
y′=f(xn,yn)，y(x0)=y0y^{\prime}= f(x_n,y_n)，y(x_0)=y_0 y′=f(xn,yn)，y(x0)=y0
则，对于该问题的RK4由如下方程给出：
yn+1=yn+h6(k1+2k2+2k3+k4)y_{n+1}=y_n+\frac h 6(k_1+2k_2+2k_3+k_4) yn+1=yn+6h(k1+2k2+2k3+k4)
其中，
k1=f(xn,yn)k_1= f(x_n,y_n) k1=f(xn,yn)
k2=f(xn+h2,yn+h2k1)k_2= f(x_n+{h\over2},y_n+{h\over2}k_1) k2=f(xn+2h,yn+2hk1)
k3=f(xn+h2,yn+h2k2)k_3= f(x_n+{h\over2},y_n+{h\over2}k_2) k3=f(xn+2h,yn+2hk2)
k4=f(xn+h,yn+hk3)k_4= f(x_n+h,y_n+hk_3) k4=f(xn+h,yn+hk3)
我开始在理解的时候花了挺长时间，毕竟公式很难让人理解，因此我才想把我遇到的困难和理解写下来，帮助大家理解龙格库塔这个方法。
在我看来龙格库塔方法就是改良之后的欧拉法，求解精度更高。就是用已知点和所求点这两个点之间几个点的斜率，加权平均后得出的斜率就为所求点与已知点的斜率。

上图就是我所举的例子，事实上我们只知晓y′y^{\prime}y′,即k1k_1k1,k2k_2k2,k3k_3k3,k4k_4k4,以及y0y_0y0。这样我们就可以通过四阶龙格库塔方法进行不断迭代，求出y1,y2…yny_1,y_2 \dots y_ny1,y2…yn。

利用四阶龙格库塔法求解一个案例

如y′=x+xyy^{\prime}=x+xyy′=x+xy这样一个表达式，我们若用正常的数学方法求解或许并不困难。但是一旦表达式变得很复杂，例如y′=sin⁡(x)cos⁡(y)+xyy^{\prime}=\sin(x)\cos(y)+xyy′=sin(x)cos(y)+xy这样的微分方程,求解起来就会相当麻烦。因此利用龙格库塔法求解，会使问题变得较简单方便。一般我们使用四阶龙格库塔方程解决问题，这样求解的精度较高，且求解速度较快。
我们举y=x+xyy=x+xyy=x+xy这样一个简单的微分方程来说明一下，由于y′=f(x,y)y^{\prime}=f(x,y)y′=f(x,y),因此f(x,y)=dydx=x+xyf(x,y)={dy \over dx}=x+xyf(x,y)=dxdy=x+xy，我们需要知道y(x0)=y0y(x_0)=y_0y(x0)=y0为多少，在此案例中,我们设x0=0x_0=0x0=0,y(0)=y0=0y(0)=y_0=0y(0)=y0=0,我们设步长h=0.1h=0.1h=0.1，计算公式如下所示：
k1=f(x0,y0)=f(0,0)=0+0=0k_1= f(x_0,y_0)=f(0,0)=0+0=0 k1=f(x0,y0)=f(0,0)=0+0=0

k2=f(x0+h2,y0+k1×h2)=f(0+0.05,0+0×0.05)=0.05+0=0.05k_2= f(x_0+\frac h 2,y_0+ k_1 \times{h \over 2})=f(0+0.05,0+0\times0.05)=0.05+0=0.05 k2=f(x0+2h,y0+k1×2h)=f(0+0.05,0+0×0.05)=0.05+0=0.05

k3=f(x0+h2,y0+k2×h2)=f(0+0.05,0+0.05×0.05)=0.05+0.05×0.0025=0.050125k_3= f(x_0+\frac h 2,y_0+k_2\times{h \over 2})=f(0+0.05,0+0.05\times0.05)=0.05+0.05\times0.0025=0.050125 k3=f(x0+2h,y0+k2×2h)=f(0+0.05,0+0.05×0.05)=0.05+0.05×0.0025=0.050125
k4=f(x0+h2,y0+k3×h2)=f(0+0.1,0+0.050125×0.1)=0.1+0.1×0.050125=0.10050125k_4= f(x_0+\frac h 2,y_0+k_3\times{h \over 2})=f(0+0.1,0+0.050125\times0.1)=0.1+0.1\times0.050125=0.10050125 k4=f(x0+2h,y0+k3×2h)=f(0+0.1,0+0.050125×0.1)=0.1+0.1×0.050125=0.10050125
y1=y0+h6(k1+2k2+2k3+k4)≈0+0.16×0.3008=0.0050y_1=y_0+\frac h 6(k_1+2k_2+2k_3+k_4)\approx0+\frac {0.1}{6} \times 0.3008= 0.0050 y1=y0+6h(k1+2k2+2k3+k4)≈0+60.1×0.3008=0.0050
到此我们已经求出y1=0.0050y_1= 0.0050y1=0.0050，就可把y1y_1y1带入下一步继续迭代，可是有的同学就会问x1x_1x1等于多少，由于我们设的步长h=0.1h=0.1h=0.1,因此x1=x0+h=0.1x_1=x_0+h=0.1x1=x0+h=0.1,扩展后则为xn=x0+nhx_n=x_0+nhxn=x0+nh。
我们需要验通过四阶龙格库塔法求解的值准不准确，因此我们用经典的微分方程解法去求解，求解过程如下：
y′=dydx=x+xyy^{\prime}=\frac {dy}{dx}=x+xy y′=dxdy=x+xy
dydx=x(1+y)\frac {dy}{dx}=x(1+y) dxdy=x(1+y)
dy1+y=xdx\frac {dy}{1+y}=xdx 1+ydy=xdx
两边同时积分
ln⁡(1+y)=12x2+C\ln{(1+y)}={1\over 2}x^2+C ln(1+y)=21x2+C
y=e12x2+Cy=e^{{1\over 2}x^2}+C y=e21x2+C
我们设x0=0x_0=0x0=0时，y0=0y_0=0y0=0,因此C=−1C=-1C=−1，即y=e12x2−1y=e^{{1\over 2}x^2}-1y=e21x2−1,当x1=0.1x_1=0.1x1=0.1时，y1=0.0050y_1=0.0050y1=0.0050。验证后我们发现四阶龙格库塔法的精确度确实很高。

用MATLAB编程

clc,clear;
syms x y %设置符号变量x,y
f(x,y)=x*y+x;
z(x)=exp(1/2*x^2)-1%我们计算得出的方程，用于验算
df=f;
y0=0;%设y0=0
a=0;%设起始点
b=3;%设终点
h=0.1;%设步长
[x1,y1]=FourLongkuta(df,y0,a,b,h);%用四阶龙格库塔法求解
figure
plot(a:h:b,z(a:h:b))
hold on
plot(a:h:b,y1,'r')
legend('原函数','龙格库塔法')
hold off

下面是四阶龙格库塔方法的函数：

function [x1,y1]=FourLongkuta(df,y0,a,b,h)%从左到右依次为已知的导函数，初始函数值，初始横坐标值，终止横坐标值，步长
n=floor((b-a)/h);%求解所需的总步长，floor的用法是取一个小数的最小整数值
y1(1)=y0;%赋值y1(1)=y(t0),且数为双精度，这样以后计算的符号值都会转换成双精度
x1=a:h:b;%设xn的值
for i=1:n%龙格库塔方法进行数值求值      k1=double (df( x1(i),y1(i) ) ) ;%double将符号值转为双精度值可以提高运算速度k2=double (df( x1(i)+h/2,y1(i)+k1*h/2 ) );k3=double (df( x1(i)+h/2,y1(i)+k2*h/2 ) );k4=double (df( x1(i)+h,y1(i)+k3*h ) );y1(i+1)=y1(i)+h*(k1+2*k2+2*k3+k4)/6;
end
end

由图可以看到四阶龙格库塔法的精度很高，几乎覆盖在原函数上。
创作不易，请大家多多支持！

龙格库塔方法的原理和案例及MTATLAB编程相关推荐

计算方法实验（三）：四阶龙格-库塔方法
四阶Runge-Kutta数学原理给定常微分方程初值问题 {dydx=f(x,y),a≤x≤by(a)=αh=b−aN\left\{ \begin{matrix} \frac{\text{dy}}{ ...
matlab:使用4阶龙格库塔方法求解常微分方程组
%书籍:常用数值算法及其matlab实现 %第10章常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
matlab 数值计算课二阶微分方程-龙格库塔方法 ODE45
详见mathworks 龙格库塔方法写成矩阵(状态方程)的形式更简洁一点(其实这两种方法结果是一样的,如果C是[1,0,0]的话,就很明显了) 例如:求系统在0-5s内的单位阶跃相应,已知传递函数: ...
matlab：使用四阶龙格库塔方法求解微分方程组
%书籍:常用数值算法及其matlab实现 %第10章常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
【Hydro】龙格-库塔方法的公式推导
一.龙格-库塔方法的公式推导摘自<数值计算方法(MATLAB版)> 二.龙格-库塔方法在水库调洪演算的应用一般库水面坡降很小,忽略动库容影响,近似看成静水面,水库蓄水量V只随坝前水位Z ...
如何用龙格库塔方法求解范德波振子的运动
如何用龙格库塔方法求解范德波振子的运动坑,没有想好. 这是一个非常好的线索: Generalized Forced Van Der Pol Oscillator Phase Plot - fun ...
用四阶龙格-库塔方法求微分方程组
最近一段时间再忙期末考试,小学期课程设计的东西,没怎么更新博客.... 更新一个用四阶龙格库塔方法求解脉冲微分方程,题目来源是一篇论文<Impulsive control of projecti ...
matlab：使用4阶龙格库塔方法求微分方程组的值
参考书目:<数值分析(matlab版)> 作者:周璐翻译 %调用龙格库塔方法求解微分方程组,P399,9.7 微分方程组 function [t,z] = my_rk4b(f, t0, ...
龙格-库塔方法学习笔记
1.龙格-库塔法简介龙格-库塔法是一种在工程上应用广泛的高精度单步算法,其中包括著名的欧拉法,用于数值求解微分方程. 由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂. 在各种 ...
二阶龙格库塔公式推导_DeepFM原理推导
注:欢迎关注本人分享有关个性化推荐系统公众号:Tiany_RecoSystem 转发一篇关于DeepFM的公式推导的博客: 论文精读-DeepFM - CSDN博客blog.csdn.net 以及D ...

龙格库塔方法的原理和案例及MTATLAB编程

文章目录

龙格库塔法的原理

利用四阶龙格库塔法求解一个案例

用MATLAB编程

龙格库塔方法的原理和案例及MTATLAB编程相关推荐

最新文章

热门文章