- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
科学计算与MATLAB
2012
第七讲 常微分方程数值解法
内容提要
引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。
例如,质点的加速运动，简谐振动等。
F m dv dt
d2x 2x2 0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
作如下近似
yn y(xn )
得：
yn1 yn f yn , xn h
2.1.4 欧拉法误差
利用泰勒级数得：
y xn1 y(xn h)
y(xn )
hy(xn )
1 2
h2
y(xn )
y(xn )
x2 …. xn ….
y(x0) y(x1) y(x2) …. y(xn) ….
y0
y1 y2 …. yn ….
在xn节点上，微分方程可以写为
y(xn1) y(xn ) f y(xn ) , xn h
作如下近似：
yn y(tn )
则得到欧拉解法递推公式的一般形式：
yn1 yn f ( yn , xn ) h
hf
y(xn ),
xn
1 2
h2 y(xn )
作如下近似
yn y(xn )
yn1 yn f yn , xn h
局部截断误差
y
y0
dy dx
x0 x x0
y0 f ( y0 , x0 ) x x0
此切线与x=x1交点纵坐标为：
y1 y0 f ( y0 , x0 ) x1 x0
=y0 f ( y0 , x0 ) h
y1
Qy((tt))
t0 t1 t2 t3 t4 t5 t6 t
从（x1,y1）作曲线y(x)在 (x1,y(x1))的切线的平行线，
作曲线y(x)在(xn,y(xn))的切
Q1
线的平行线，得直线方程：
Qy((tt))
y yn f ( yn , xn ) x xn
与x=xn+1交点纵坐标为：
tx0 1 tx1 2 tx23 xt34 tx45 tx56 t6x7 t
yn1 yn f ( yn , xn ) h
折线近似曲线，n增大，误差变大
Qy
得x1
y1 f ( y1, x1) x x1
与x=x2交点纵坐标为：
y2 y1 f ( y1, x1) h
y1
Qy((tt))
y2
tx0 1 tx1 2 tx23 xt34 tx45 tx56 t6x7 t
y1 y(x1)
按照相似的方法，从（xn,yn）Qy
yn y(xn )
2、欧拉近似方法
2.1 简单欧拉(L.Euler, 1707-1783)方法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一个函数值由前一个函数值得到。关 键是构造递推公式。
y0 y1 y2 yn
欧拉数值算法递推公式构造 2.1.1 差分法
f=f(:)'; y(k + 1,:) =y(k,:) +h*f; %对于所取的点x迭代计算y值
end
outy=y;
outx=x; %plot(x,y)%画出方程解的函数图
例题： y ' sin x y y(x0 ) 1, x0 0
2.1.2 折线法（几何意义） Qy 从（x0,y0）作曲线y(x)的切 线，得切线方程：
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)在每个迭代点的值
同样，在[x0,x2] ，积分采用矩形近似，得：
y(x2 ) y0
x2 f y(x), x dx
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
2.1.3 数值积分
对微分方程作从 x0 到 x 积分得：
y(x) dy y0
x
x0 f y
x
, xdx
y x y0
x x0
f y x, xdx
在[x0,x1] ，积分采用矩形近似，得：
y(x1) y0
x1 f y(x), x dx
x0
y0 f y(x0 ), x0 h
dt 2
简单问题可以求得解析解，多数实际问题靠数值求解。
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立的节点 xn 上的近似值 yn，用以近似y(xn)
差分法就是用差商近似代替微商，即
y dy x dx
代入微分方程得到：
y(x x) y(x) f ( y(x), x) x
y(x x) y(x) f ( y(x), x)x
对于等间隔节点
x xn1 xn h xn1 xn h
n=0,1,2
可以得到：
xn y精确值 y近似值
x0 x1
具体求解过程为：
y1 y0 f ( y0 , x0 ) h y2 y1 f ( y1 , x1) h y3 y2 f ( y2 , x2 ) h
简单欧拉方法程序
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
