常微分方程初值问题的Runge-Kutta方法.docVIP

1011203002 自动化学院 李霞林 PAGE 3 常微分方程初值问题的Runge-Kutta方法 背景 在自然科学和工程技术中的许多问题都可以用微分方程描述。设常微分方程（或称动力系统） （1） 大量工程技术问题中所出现的微分方程往往不能直接求出解析解，因此采用数值方法求解近似解，使得（1）积分轨道上tn+1时刻一个点的xn+1可以通过tn或前几个时刻的x值（也称步点）来获得。若xn+1通过前面一个步点求得，这样的数值方法称单步法，单步法中最常用的是Runge-Kutta方法。 算法思想 考虑一阶常微分方程初值问题 （2） 其中，f ( tn , un ) 是区域 D = 上的实值函数。 Taylor展开式表明用函数在一个点上的各阶导数可以近似表示它在邻近一个点上的函数值；数值微分公式表明函数在一个点上的各阶导数值可以用邻近一些点上的函数值近似表示。RK方法基于后一种思路，用f(t,u)在一些特殊点上函数值的线性组合来表示Taylor展开式的各阶导数值。经典的四阶Runge-Kutta算法为 (3) 实例与编程 利用经典4阶Runge-Kutta格式求解初值问题 （4） 由式 (3) ，可以写出上述一阶微分方程的标准Runge-Kutta格式为 （5） 取步长为 h = 0.4 ， Matlab程序为 clear; h=0.4; t(1)=0; u(1)=1; for i=1:2/h k1=u(i)+2*t(i)/u(i)^2; k2=u(i)+0.5*h*k1+2*(t(i)+0.5*h)/(u(i)+0.5*h*k1)^2; k3=u(i)+0.5*h*k2+2*(t(i)+0.5*h)/(u(i)+0.5*h*k2)^2; k4=u(i)+1.0*h*k3+2*(t(i)+1.0*h)/(u(i)+1.0*h*k3)^2; t(i+1)=h*i; u(i+1)=u(i)+(h/6)*(k1+2*k2+2*k3+k4); end 结果分析 表1给出了Euler法和4阶RK法对算例的数值结果，从对比结果可以看出两 者差别较大。 表1 计算结果 tn Euler un 四阶RK法 un 0 1.0000 1.0000 0.4 1.4000 1.5963 0.8 2.1233 2.5254 1.2 3.1145 3.8693 1.6 4.4593 5.8346 2.0 6.3074 8.7395 利用matlab自带函数dslove()可以求解出上述算例的精确解，求解方式如下： syms x y y=dsolve(Dy-y-2*x/y^2=0,y(0)=1,x) ezplot(y) axis([0 2 1 9]) 可以求解得到式（4）的精确解为： （6） 利用RK法求出的结果画图，如图1所示，利用精确解得到的结果如图2所示。可以看出RK法求得的数值解精度要比Euler的高。 图1 图2