matlab_常微分方程数值解法

简单问题可以求得解析解，多数实际问题靠数值求解。
2
一阶常微分方程(ODE )初值问题 ：
ODE ：Ordinary Differential Equation
dy f ( x, y ) x0 x xn dx y ( x0 ) y0
数值解法就是求y(x)在某些分立的节点 xn 上的近似值 yn，用以近似y(xn)
在[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 向后的欧拉方法 取后一点的斜率作为平均斜率，即：
yn1 yn k2 h k2 f ( xn1 , yn1 )
整体截断误差： 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
(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 )
对于等间隔节点
x xn1 xn h xn1 xn h
可以得到： xn y精确值 y近似值 x0 x1 x2
n=0,1,2
…. …. ….
xn
…. ….
y(x0) y(x1) y(x2) y0 y1 y2
y(xn) yn
….
在xn节点上，微分方程可以写为
y( xn1 ) y( xn ) f y( xn ) , x n h
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( xn )
2、欧拉近似方法
2.1 简单欧拉(L.Euler, 1707-1783)方法。
dy f ( y, x) dx y ( x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解 就是从初值开始，后一个函数值由前一个函数值得到。关 键是构造递推公式。


改进的欧拉方法递推公式为
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
按照相似的方法，从（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
作如下近似
yn y( xn )
yn1 yn f yn , xn h
2.2 改进的欧拉近似方法 对于一般的常微分方程：
dy f ( x, y ) dx y ( x0 ) y0
节点： 步长：
x0 , x1 ,, xi ,, xn , h xi 1 xi
格式 [x,y]=ode45(′ fun′, [x0,xn], y0,option]
说明： yi( x) fi ( x, yi ) ① 适用于求解一阶常微分方程组 ② fun定义微分方程组的函数文件名 ③ [x0,xn]求解区域 ④ y0初始条件向量 ⑤ option可选参数，由ODESET函数设置，比较复杂 ⑥ x输出自变量向量，y输出[y, y ′, y ″,..]
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
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
)
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
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
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 xn 1 , y xn 1 2
例题：
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
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
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
科学计算与MATLAB
常微分方程数值解法
内容提要



引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。 例如,质点的加速运动，简谐振动等。
dv F m dt
d x 2 2 x 0 2 dt
整体截断误差： O(h4)
例题： y ' sin x y
y ( x0 ) 1, x0 0
time_Euler = 0.2500 time_EulerPro = 0.5000 time_RK4 = 1.0150
4、 MATLAB的常微分方程函数
ode45 ode23 ode113 ode15s ode23s ode23t ode23tb
简单欧拉方法程序
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)%画出方程解的函数图