4阶RungeKutta法求解一阶常微分方程Word文档格式.docx

4阶RungeKutta法求解一阶常微分方程Word文档格式.docx

《4阶RungeKutta法求解一阶常微分方程Word文档格式.docx》由会员分享，可在线阅读，更多相关《4阶RungeKutta法求解一阶常微分方程Word文档格式.docx（6页珍藏版）》请在冰豆网上搜索。

龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。

如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。

一阶常微分方程可以写作：

y'

=f（x,y）,使用差分概念。

（Yn+1-Yn）/h=f（Xn,Yn）推出（近似等于，极限为Yn'

）Yn+1=Yn+h*f（Xn,Yn）

另外根据微分中值定理，存在0<

t<

1,使得Yn+1=Yn+h*f（Xn+th,Y（Xn+th））

这里K＝f（Xn+th,Y（Xn+th））称为平均斜率，龙格库塔方法就是求得K的一种算法。

利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O（h^5）的四阶龙格库塔公式：

K1＝f（Xn,Yn）;

K2=f（Xn+h/2,Yn+（h/2）*K1）;

K3=f（Xn+h/2,Yn+（h/2）*K2）;

K4=f（Xn+h,Yn+h*K3）;

Yn+1=Yn+h*（K1+2K2+2K3+K4）*（1/6）

二、Runge-Kutta的算法和流程图

在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为0（h5）,被广泛应用于解微分方程的初值问题。

其算法公式为:

流程图：

（1）、四阶龙格-库塔方法流程图:

（2）、实例求解流程图：

三、Runge-Kutta的Matlab实现

function[x,y]=runge_kutta1（ufunc,y0,h,a,b）%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor（（b-a）/h）;

%求步数x

（1）=a;

%时间起点y（:

1）=y0;

%赋初值，可以是向量，但是要注意维数forii=1:

n

x（ii+1）=x（ii）+h;

k1=ufunc（x（ii）,y（:

ii））;

k2=ufunc（x（ii）+h/2,y（:

ii）+h*k1/2）;

k3=ufunc（x（ii）+h/2,y（:

ii）+h*k2/2）;

k4=ufunc（x（ii）+h,y（:

ii）+h*k3）;

y（:

ii+1）=y（:

ii）+h*（k1+2*k2+2*k3+k4）/6;

%按照龙格库塔方法进行数值求解end

四、Runge-Kutta的算例实现

例:

求解常微分方程d（y）/d（x）=-2*y+2*x*x+2*x，0≤x≤0.5，y（0）=1.

编辑m文件Untitled3：

fun=inline（'

-2*y+2*x*x+2*x'

）;

[x,y]=ode45（fun,[0,0.5],1）

>

Untitled3

x=

0

0.0125

0.0250

0.0375

0.0500

0.0625

0.0750

0.0875

0.1000

0.1125

0.1250

0.1375

0.1500

0.1625

0.1750

0.1875

0.2000

0.2125

0.2250

0.2375

0.2500

0.2625

0.2750

0.2875

0.3000

0.3125

0.3250

0.3375

0.3500

0.3625

0.3750

0.3875

0.4000

0.4125

0.4250

0.4375

0.4500

0.4625

0.4750

0.4875

0.5000

y=

1.0000

0.9755

0.9519

0.9291

0.9073

0.8864

0.8663

0.8471

0.8287

0.8112

0.7944

0.7785

0.7633

0.7489

0.7353

0.7224

0.7103

0.6989

0.6883

0.6783

0.6690

0.6605

0.6526

0.6454

0.6388

0.6329

0.6277

0.6231

0.6191

0.6157

0.6130

0.6109

0.6093

0.6084

0.6080

0.6083

0.6091

0.6104

0.6124

0.6148

0.6179