常微分方程

Ordinary differential equation，简称ODE，自变量只有一个的微分方程。
例子1： dydx=f(x,y)\dfrac {dy} {dx}=f(x,y)dxdy=f(x,y) ,f(x,y)f(x,y)f(x,y)是已知函数

偏微分方程

Partial differential equation，简称PDE，自变量有多个的微分方程。
例子2：ut−a2uxx=0,a>0u_t-a^2u_{xx}=0,a>0ut−a2uxx=0,a>0为常数（热传导方程，抛物型方程的典型代表）

显式（Explicit）

第n步结果可以从n-1, n-2, …1步的结果直接推导出来，迭代时每步的计算量很小，但迭代增量也有限制，不能太大，否则会出现发散。

隐式（Implicit）

第n步的计算结果不能直接从前面的结果推导出来，必须做进一步的求解，这样，迭代时每步的计算量很大。

泰勒公式

f(x)=f(a)0!+f′(a)1!(x−a)+f′′(a)2!(x−a)2+...+f(n)(a)n!(x−a)n+Rn(x)f(x)=\dfrac {f(a)} {0!}+\dfrac {f'(a)} {1!}(x-a)+\dfrac {f''(a)} {2!}(x-a)^2+...+\dfrac {f^{(n)}(a)} {n!}(x-a)^n+R_n(x) f(x)=0!f(a)+1!f′(a)(x−a)+2!f′′(a)(x−a)2+...+n!f(n)(a)(x−a)n+Rn(x)
称为 f(x) 在 x = a 点关于 x 的幂函数展开式，又称为 Taylor 公式，式中Rn（x）叫做 Lagrange 余项。

欧拉方法

考虑一阶常微分方程的初值问题
{dydx=f(x,y)y(x0)=y0\left\{ \begin{aligned} &\dfrac {dy} {dx}=f(x,y)\\ &y(x_0)=y_0 \end{aligned} \right.⎩⎨⎧dxdy=f(x,y)y(x0)=y0

前向欧拉法

yn+1=yn+hf(xn,yn),n=0,1,...y_{n+1}=y_n+hf(x_n,y_n),n=0,1,...yn+1=yn+hf(xn,yn),n=0,1,...

后向欧拉算法

yn+1=yn+hf(xn+1,yn+1),n=0,1,...y_{n+1}=y_n+hf(x_{n+1},y_{n+1}),n=0,1,...yn+1=yn+hf(xn+1,yn+1),n=0,1,...
前后向欧拉法的推理、举例及稳定性对比的超棒分析！地址入口

梯形公式

yn+1=yn+h2[f(xn+yn)+f(xn+1+yn+1)],n=0,1,...y_{n+1}=y_n+\frac h 2[f(x_n+y_n)+f(x_{n+1}+y_{n+1})],n=0,1,...yn+1=yn+2h[f(xn+yn)+f(xn+1+yn+1)],n=0,1,...

改进的欧拉算法

{y^n+1=yn+hf(xn,yn)yn+1=yn+h2[f(xn,yn)+f(xn+1,y^n+1)]\left\{ \begin{aligned} &\hat{y}_{n+1}=y_n+hf(x_n,y_n)\\ &y_{n+1}=y_n+\frac h 2[f(x_n,y_n)+f(x_{n+1},\hat{y}_{n+1})] \end{aligned} \right.⎩⎨⎧y^n+1=yn+hf(xn,yn)yn+1=yn+2h[f(xn,yn)+f(xn+1,y^n+1)]
也可以表示成
{yf=yn+hf(xn,yn)yb=yn+hf(xn+1,yf)yn+1=12(xf+xb)\left\{ \begin{aligned} &y_f=y_n+hf(x_n,y_n)\\ &y_b=y_n+hf(x_{n+1},y_f)\\ &y_{n+1}=\frac 1 2 (x_f+x_b) \end{aligned} \right.⎩⎪⎪⎪⎨⎪⎪⎪⎧yf=yn+hf(xn,yn)yb=yn+hf(xn+1,yf)yn+1=21(xf+xb)

其中yfy_fyf表示利用向前(显式)欧拉公式的近似值，yby_byb表示利用向后(隐式)欧拉公式的近似值(利用了yfy_fyf)，最后取平均值。
上面利用f(xn+1,y^n+1)]f(x_{n+1},\hat{y}_{n+1})]f(xn+1,y^n+1)]是有误差的，也可通过多次迭代减少误差，具体公式如下
{y^n+1(0)=yn+hf(xn,yn)yn+1(k+1)=yn+h2[f(xn,yn)+f(xn+1,y^n+1(k))],k=0,1,2,...\left\{ \begin{aligned} &\hat{y}_{n+1}^{(0)}=y_n+hf(x_n,y_n)\\ &y_{n+1}^{(k+1)}=y_n+\frac h 2[f(x_n,y_n)+f(x_{n+1},\hat{y}_{n+1}^{(k)})],k=0,1,2,... \end{aligned} \right.⎩⎪⎨⎪⎧y^n+1(0)=yn+hf(xn,yn)yn+1(k+1)=yn+2h[f(xn,yn)+f(xn+1,y^n+1(k))],k=0,1,2,...

其中，后向欧拉算法和梯形公式是隐式算法，前向欧拉算法和改进的欧拉算法是显式算法。欧拉算法计算容易，但是精度低，梯形公式精度高，但是是隐式形式，不易求解。将两式结合，则可以得到改进的欧拉公式。先用欧拉公式求出yn+1的一个粗糙的估计值，再用梯形方法进行精确化，称为校正值。

四阶龙格库塔算法（Runge-kutta method）

{k1=f(xn,yn)k2=f(xn+h2,yn+h2×k1)k3=f(xn+h2,yn+h2×k2)k4=f(xn+h,yn+h×k3)yn+1=yn+(k1+2k2+2k3+k4)6×h\left\{ \begin{aligned} &k_1=f(x_n,y_n)\\ &k_2=f(x_n+\frac h 2,y_n+\frac h 2×k_1)\\ &k_3=f(x_n+\frac h 2,y_n+\frac h 2×k_2)\\ &k_4=f(x_n+h,y_n+h×k_3)\\ &y_{n+1}=y_n+\frac {(k_1+2k_2+2k_3+k_4)} 6 ×h \end{aligned} \right.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧k1=f(xn,yn)k2=f(xn+2h,yn+2h×k1)k3=f(xn+2h,yn+2h×k2)k4=f(xn+h,yn+h×k3)yn+1=yn+6(k1+2k2+2k3+k4)×h
龙格-库塔算法是一种在工程上应用广泛的高精度单步算法，算法精度较高，如果预先取四个点就是四阶龙格-库塔算法。
真的！龙格库塔算法的超强解析！地址入口

四阶龙格库塔函数代码

runge_kuttx0_o4.m
代码参考的博客地址入口
总结一下代码的优点~
1、使用函数句柄，方便复用~
2、模仿了ode45的函数输入变量，和ode45用起来差不多~

function [t,y,n]=runge_kuttx0_o4(ufunc,tspan,y0,h)%参数表顺序依次是微分方程组的函数名称，时间起点终点,初始值，步长（参数形式参考了ode45函数）
if nargin<4h=0.01;
end
if size(tspan)==[1,2]t0=tspan(1);tn=tspan(2);
elseerror(message('MATLAB:runge_kuttx0_o4:WrongDimensionOfTspan'));
end
n=floor((tn-t0)/h);%求步数
t(1)=t0;%时间起点
y(:,1)=y0;%第一列赋初值，可以是向量，代表不同的初始值
for i=1:n
t(i+1)=t(i)+h;
k1=ufunc(t(i),y(:,i));
k2=ufunc(t(i)+h/2,y(:,i)+h*k1/2);
k3=ufunc(t(i)+h/2,y(:,i)+h*k2/2);
k4=ufunc(t(i)+h,y(:,i)+h*k3);
y(:,i+1)=y(:,i)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end

test1.m

clear
clc
test_fun=@(t,y)(y+3*t)/t^2;
tspan=[1 4];
y0=-2;
h=1;
[t1,y1]=ode45(test_fun,tspan,y0);
[t2,y2]=runge_kuttx0_o4(test_fun,tspan,y0,h);
plot(t1,y1,'r',t2,y2,'g')
legend('ode45函数效果','自编四阶龙格库塔函数效果')
xlabel('t');
ylabel('y');
title('效果对比图')

test.m运行结果

步长选择

之前讨论的所有的龙格-库塔方法都是以ΔtΔtΔt定步长来展开的，但从xi⇒xi+1x_i ⇒x_{i+1}xi⇒xi+1单步递推过程来说，步长ΔtΔtΔt越小，局部截断误差越小(方法确定情况下)，但是随着步长的缩小，不但会引起计算量的增加，而且也有可能引起舍入误差的严重积累；但步长ΔtΔtΔt太大又不能达到预期的精度要求，所以选择合适的步长ΔtΔtΔt，在实际计算中也是比较重要的。其实有时候在实际使用中步长并不需要算法确定，而是需要根据数据帧率来确定的，比如imu数据。??
下面给出求解步长的步骤：

1、以步长ΔtΔtΔt开始，利用龙格-库塔公式计算xi⇒xi+1x_i ⇒x_{i+1}xi⇒xi+1得到一个近似值xi+1Δtx_{i+1}^{\Delta t}xi+1Δt；
2、然后步长减半为Δt/2\Delta t /2Δt/2，利用龙格-库塔公式分两步计算xi⇒xi+12⇒xi+1x_i ⇒ x_{i + \frac1 2} ⇒ x_{i+1}xi⇒xi+21⇒xi+1得到一个近似值xi+1Δt/2x_{i + 1}^{Δt/2}xi+1Δt/2；
3、计算∣xi+1Δt/2−xi+1Δt<ϵ∣∣x_{i + 1}^{Δt/2}-x_{i + 1}^{Δt}< ϵ ∣∣xi+1Δt/2−xi+1Δt<ϵ∣是否成立，如果成立直接步长选择ΔtΔtΔt，否则继续步长减半，重复上诉步骤直到满足精度要求。
test2.m

clear
clc
test_fun=@(t,y)(y+3*t)/t^2;
tspan=[1 15];
y0=-2;
h=1;
acc=0.0001;
[t1,y1,n1]=runge_kuttx0_o4(test_fun,tspan,y0,h);
tf=t1;yf=y1;hf=h;
h=h/2;
[t2,y2,n2]=runge_kuttx0_o4(test_fun,tspan,y0,h);
while(sum(abs(y1(2:n1+1)-y2(3:2:n2+1)))/n1>=acc)
t1=t2;y1=y2;n1=n2;
h=h/2;
[t2,y2,n2]=runge_kuttx0_o4(test_fun,tspan,y0,h);
end
plot(tf,yf,'r',t2,y2,'g')
legend(['四阶龙格库塔初始步长下的效果,h=',num2str(hf)],['四阶龙格库塔精度约为0.01的效果,h=',num2str(h*2)])
xlabel('t');
ylabel('y');
title('效果对比图')

四阶龙格库塔算法及matlab代码相关推荐

四阶龙格库塔算法用MATLAB写
四阶龙格-库塔算法可以使用 MATLAB 进行编写.您可以使用 MATLAB 的数值解法工具箱来解决常微分方程组,并使用相应的函数(例如 ode45)来实现四阶龙格-库塔算法.在编写代码时,您需要根据 ...
四阶龙格库塔c语言,四阶龙格库塔算法的C语言实现
解微分方程 2001年3月焦作大学学报 JOURNALOFJIAOZUOUNIVERSITY№ 1Mar．2001第1期四阶龙格一库塔算法的C语言实现毋玉芝 (焦作财会学校) 摘要本文叙述了四阶龙 ...
四阶龙格库塔方程解二阶常微分方程组并计算船舶在迎浪下的纵摇埀荡耦合运动方程-附Python代码
0 写在前面这篇博客是在将我上一篇博客的matlab代码移植到python中,应为后续要开展深度强化学习下的船舶减摇研究,总的来说还是在python上做这项工作比较合适.具体方程的解法和背景不在赘述 ...
捷联惯导基础知识之一（姿态更新的毕卡算法、龙格库塔算法、及精确数值解法）
一.介绍旋转矢量以及由旋转矢量更新四元数等效旋转矢量微分方程(Bortz方程): 经过理论和近似推导: 单位四元数与旋转矢量关系: 单位四元数与三角函数的关系: 四元数更新方程:利用旋转矢量,考虑了 ...
四阶龙格库塔方程（Rungekutta）解二阶常微分方程组并计算船舶在迎浪下的纵摇埀荡耦合运动方程-附Matlab代码
今年年初的时候给师姐做了DDPG算法的船舶减横摇控制算法,师姐还有想法要让我把纵摇-埀荡两个自由度的减摇也做出来,这个任务归我了.实际上不管是多少个自由度的减摇,其实都需要解运动方程,当初做单自由度横 ...
一阶欧拉近似matlab,MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程.doc
MATLAB改进欧拉法与四阶龙格-库塔求解一阶常微分方程姓名:樊元君学号:2012200902 日期:2012.11.06 一.实验目的掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题 ...
matlab：使用四阶龙格库塔方法求解微分方程组
%书籍:常用数值算法及其matlab实现 %第10章常微分方程初值问题的数值解法,例10.14使用 %四阶龙格库塔方法 function [t,z] = rk4symeq(fun, t0, tf, ...
二阶水箱 matlab 四阶龙格库塔,请问这个二阶常微分方程组用龙格四阶库塔法怎么编写...
要是不用MATLAB自带的ode45函数也可以网上下载一个4阶龙格库塔算法来代替附上一个仅供参考 function y=DELGKT4_rungekuta(f,h,a,b,y0,varvec) ...
Ode45以及龙格-库塔算法
Ode45及龙格-库塔算法 ODE45 ODE45函数描述: 求解非刚性微分方程 - 中阶方 ODE是Matlab专门用于解微分方程的功能函数.该求解器有变步长(variable-step)和定步长( ...
数值分析C++实现用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题
问题:用四阶龙格-库塔(Runge-Kutta)方法求解常微分方程初值问题. 算法描述算法的程序框图: 具体算法: (1)读取a,b,n,f (2)计算步长h = (b-a)/n,x0=a,y0=f ...

四阶龙格库塔算法及matlab代码