实验Euler法求微分方程的解.ppt

数学实验
实验四
求微分方程的解
问题背景和实验目的
口自牛顿发明微积分以来,微分方程在描述事物运
动规律上已发挥了重要的作用。实际应用问题通过
数学建模所得到的方程,绝大多数是微分方程。
口由于实际应用的需要,人们必须求解微分方程。
然而能够求得解析解的微分方程十分有限,绝大多
数微分方程需要利用数值方法来近似求解。
口本实验主要研究如何用 Matlab来计算微分方程
(组)的数值解,并重点介绍一个求解微分方程的
基本数值解法一一 Euler折线法
Euler折线法
口考虑一维经典初值问题
f(r, y), y(o=yo, Ela,b]
d
◆基本思想:用差商代替微商
根据 Talyor公式,y(x)在点xk处有
y(r)=y(xk)+(x-xry'(xk)+O(Ar)
Ay(+)=y(xx)+hy'()+0(h2) h=xk+l-Xk
y(rk-y(xk)
+O(h)≈
dr
h
k
初值问题的Euer折线法
口具体步骤:分割求解区间,差商代替微商,解代数方程
◆分割求解区间
x}%为分割点
等距剖分:a=x0<x1<x
<∴<x.,<x.=b
步长:h=xk
xk=(b-a)/n,k=O,1,2,,n
◆差商代替微商
y(r
h
口yx+)=y(x)+hy(x)
yo=y(xo)
k=0,1,2,…,n
得方程组y=y+h(x,y)
h
是y(xk)的近似
Euler折线法举例
例:用 Euler法解初值问题
d r+<t
x∈[O,2]
y(O)=1
解:取步长h=(2-0)/n=2/m,得差分方程
yu=k+hf(rk, yk)=yk+h(yk +2xk/ yk)
x=xth
当h=04,即n=5时, Matlab源程序见
Euler折线法源程序
clear
f=sym(y+2*x/y^21);
=0;b=2;

n=(b-a)/h+1;8n=(b-a)/h;
x=0:y=1
szj=[XY]
for i=1: n-11=1: n
y=y+h*subs(f,{'x!,"'y'},{x,y})
xx+h
sz]=[szily]
sZJ
p1ot(szj(:;1),szj(:2),"or-1)
Euler折线法举例(续)
解析解:y=1/3*(-18-54*x+45*(3*x)~(1/3)
解析解
近似解
Runge-Kutta法
口为了减小误差,可采用以下方法:
◆让步长h取得更小一些
◆改用具有较高精度的数值方法:
Runge- Kutta(龙格-库塔)方法
口龙格-库塔方法
◆是一类求解常微分方程的数值方法
◆有多种不同的迭代格式
Runge-Kutta法
口用得较多的是四阶RK方法(教材第98页
yo =y(r
k+1=x+h
yk+=yk+h(1+2L2+2L3+L46
其中
L=f(rk, yr)
L,2=f(x+h/,yx +hL, /2)
L =f(x+h/2, yx +hL,/2)
L 4=f(xk+,yx+hl3)
四阶RK方法源程序
clear
f=sym(y+2*x/y^21);
a=0;b=2;h=;
n=(b-a)/h+1;n=(b-a)/h;
SZJ=[xY]
fi=1:n-18i=1:n
11=subs(f,{'x',"y'},{x,y})
12=subs(￡,{'x',y'},(x+h/2,y+11*h/2});
13=subs(f,{'x',"y'},(x+h/2,y+12*h/2});
14=subs(f,{'x!,"y'},{x+h,y+13*h})
y=y+h*(11+2*12+2*13+14)/6
xex+h
SZJ=[sz]ix,]
end
p1ot(szj(:,1),sz(:,2),'dg