matlab 常微分方程数值解法

简单欧拉方法程序
function [outx,outy]=MyEuler(fun,x0,xt,y0,PointNum) %MyEuler 用前向差分的欧拉方法解微分方程 %fun 表示f（x，y） %x0,xt表示自变量的初值和终值 %y0表示函数在x0处的值,其可以为向量形式 %PointNum表示自变量在[x0,xt]上取的点数 if nargin<5 | PointNum<=0 %如果函数仅输入4个参数值，则PointNum默认值为100 PointNum=100; end if nargin<4 %y0默认值为0 y0=0; end h=(xt-x0)/PointNum;%计算步长h x=x0+[0:PointNum]'*h;%自变量数组 y(1,:) = y0(:)';%将输入存为行向量，输入为列向量形式 for k = 1:PointNum f=feval(fun,x(k),y(k,:));%计算f(x,y)在每个迭代点的值 f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值 end outy=y; outx=x; %plot(x,y)%画出方程解的函数图
数值解法就是求y(x)在某些分立的节点 xn 上的近似值 yn，用以近似y(xn)
yn y( xn )
2、欧拉近似方法
2.1 简单欧拉(L.Euler, 1707-1783)方法。
dy f ( y, x) dx y ( x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解 就是从初值开始，后一个函数值由前一个函数值得到。关 键是构造递推公式。
在[x0,x1] ，积分采用矩形近似，得：
y ( x1 ) y0 f y ( x), x dx
x1
y0 f y ( x0 ), x0 h
x0
同样，在[x0,x2] ，积分采用矩形近似，得：
y ( x2 ) y0 f y ( x ), x dx
整体截断误差： O(h)
2 向后的欧拉方法 取后一点的斜率作为平均斜率，即：
yn 1 yn k2 h k2 f ( xn 1 , yn 1 )
整体截断误差： O(h)
3 改进的欧拉方法 取两点斜率平均值，即：
h y y k k n 1 n 1 2 2 k1 f ( xn , yn ) k f ( x , y ) n 1 n 1 2
科学计算与MATLAB
主讲：唐建国 中南大学材料科学与工程学院 2012
第七讲
常微分方程数值解法
内容提要



引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。 例如,质点的加速运动，简谐振动等。
x2 x0
y0 f y ( x), x dx f y ( x), x dx
x1 x2
y ( x1 ) f y ( x1 ), x1 h
x0
x1
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y ( xn 1 ) y0
对于等间隔节点
x xn1 xn h xn1 xn h
可以得到： xn y精确值 x0 x1 x2
n=0,1,2
…. …. ….
xn
…. ….
y(x0) y(x1) y(x2)
y(xn)
y近似值
y0
y1
y2
yn
….
在xn节点上，微分方程可以写为
y( xn1 ) y( xn ) f y( xn ) , x n h
dv F m dt
d x 2 2 x 0 2 dt
简单问题可以求得解析解，多数实际问题靠数值求解。
2
一阶常微分方程(ODE )初值问题 ：
ODE ：Ordinary Differential Equation
dy f ( x, y ) x0 x xn dx y ( x0 ) y0
即

y y0
dy f ( x, y)dx
x0 xn1 xn
x
y( x) y0 f ( x, y)dx
x0
x
y ( xn 1 ) y ( xn )
f ( x, y )dx
y( xn ) h f xn 1 , y xn 1
向后的欧拉方法递推公式为
y1
Q(t) y(t)
y2
dy y y1 x1 x x1 dx y1 f ( y1 , x1 ) x x1
与x=x2交点纵坐标为：
tx t34 tx x x x7 t 0 1 tx 12 t 56 t6 23 x 45 t
y1 y( x1 )
y2 y1 f ( y1 , x1 ) h
按照相似的方法，从（xn,yn） Q y 作曲线y(x)在(xn,y(xn))的切 线的平行线，得直线方程：
Q1
y(t) Q(t)
y
yn f ( yn , xn ) x xn
与x=xn+1交点纵坐标为：
tx t34 tx x x x7 t 0 1 tx 12 t 56 t6 23 x 45 t
2.2 改进的欧拉近似方法 对于一般的常微分方程：
dy f ( x, y ) dx y ( x0 ) y0
节点： 步长：
x0 , x1 ,
, xi ,
, xn ,
h xi 1 xi
积分法求欧拉公式时，矩形法采用前一点函数值求积， 若利用后一点，即为向后的欧拉方法。
整体截断误差： O(h4)
例题： y ' sin x y
y ( x0 ) 1, x0 0
time_Euler = 0.2500 time_EulerPro = 0.5000 time_RK4 = 1.0150
t0 t1 t2
t3 t4 t5 t6
t
y1 y0 f ( y0 , x0 ) x1 x0 =y0 f ( y0 , x0 ) h
此切线与x=x1交点纵坐标为：
y0 y( x0 )
从（x1,y1）作曲线y(x)在 Q y (x1,y(x1))的切线的平行线， 得直线方程：
作如下近似
yn y( xnBaidu Nhomakorabea)
yn1 yn f yn , xn h
局部截断误差 由泰勒展式可以看出，在欧拉方法中，假定yn 是精 确的，由yn计算yn+1所产生的误差的数量级为O(h2) 。 整体截断误差 如果再考虑到yn本身的误差， 计及误差的累积效果 和局部截断误差，这种误差称为欧拉方法整体截断误差， 数量级为O(h) 。 欧拉法是数值方法解微分方程最简单的一种，精度 不高，实际应用不多，但它可以反映数值求解微分方程 的特征。
yn1 yn f ( yn , xn ) h
折线近似曲线，n增大，误差变大
2.1.3 数值积分
对微分方程作从 x0 到 x 积分得：

y( x ) y0
dy f y x , x dx x0
x
x
y x y0 f y x , x dx x0
(0) n 1
例题： y ' sin x y
y ( x0 ) 1, x0 0
3、龙格-库塔(R-K)方法
几种方法的比较 对于一般的常微分方程：
dy f ( x, y ) dx y ( x0 ) y0
1 欧拉方法 取前一点的斜率作为平均斜率，即：
yn 1 yn k1 h k1 f ( xn , yn )
y0 y1 y2
yn
欧拉数值算法递推公式构造
2.1.1 差分法 差分法就是用差商近似代替微商，即
y dy x dx
代入微分方程得到：
y ( x x) y ( x) f ( y ( x), x) x y ( x x) y ( x) f ( y ( x), x)x
yn y( xn )
yn1 yn h f ( xn1 , yn1 )
向后的欧拉方法(隐式方法):预报---校正法
1. 用欧拉方法预报 2. 用向后的欧拉方法校正 ① 用欧拉方法预报
y
(0) n1
yn h f ( xn , yn )
② 用向后的欧拉方法校正---迭代
xn1
y ( xn ) f y ( xn ), xn h
作如下近似
x0
f y ( x), x dx
yn y( xn )
得：
yn1 yn f yn , xn h
2.1.4 欧拉法误差
利用泰勒级数得：
y xn 1 y ( xn h) 1 2 y ( xn ) hy( xn ) h y( xn ) 2 1 2 y ( xn ) hf y ( xn ), xn h y( xn ) 2
作如下近似：
yn y(tn )
则得到欧拉解法递推公式的一般形式：
yn1 yn f ( yn , x n ) h
具体求解过程为：
y1 y0 f ( y0 , x0 ) h y2 y1 f ( y1 , x1 ) h
y3 y2 f ( y2 , x 2 ) h
例题：
y ' sin x y y ( x0 ) 1, x0 0
2.1.2 折线法（几何意义）
Q y
从（x0,y0）作曲线y(x)的切 线，得切线方程：
y1
y(t) Q(t)
dy y y0 x0 x x0 dx y0 f ( y0 , x0 ) x x0


改进的欧拉方法递推公式为
h yn 1 yn f ( xn , yn ) f ( xn 1 , yn 1 ) 2
隐式方法: 预报---校正法
1. 用欧拉方法预报
2. 用改进的欧拉方法校正
计算公式为：
y yn h f ( xn , yn ) h (0) f ( x , y ) f ( x , y ) yn 1 yn n n n 1 n 1 2
y
( k 1) n1
yn h f ( xn1, y )
(k ) n1
2.2.2 改进的欧拉方法 积分法求欧拉公式时，积分采用梯形近似，即得到改 进的欧拉方法。
y ( xn 1 ) y ( xn )
xn1
xn
f ( x, y )dx
h y ( xn ) f xn , y xn f xn1 , y xn 1 2
整体截断误差： O(h2)
欧拉方法（前、后）和改进的欧拉方法优点是算法简 单，但计算精度较低，不能满足实际问题求解对精度的要 求，所以使用较少。 龙格-库塔(R-K)方法就是一种精度较高的较为实用计算 常微分方程的方法。
从以上三种方法的比较，可以推测：若取多点处斜 率的加权平均作为平均斜率，精度会更高，这就是龙 格—库塔法的基本思想。
yn 1 yn ci ki
i 1
N
k1 f ( xn , yn ) i 1 ki f ( xn ai h, yn bij k j ) j 1
四阶龙格—库塔法计算公式为：
h yn 1 yn k1 2k2 2k3 k4 6 k1 f ( xn , yn ) h h k2 f ( xn , yn k1 ) 2 2 h h k f ( x , y k ) 3 n n 2 2 2 k f ( x h, y h k ) n n 3 4