第六章 常微分方程初值问题初步1

dy
function E=euler(f,a,b,ya,N) %f为该问题中的函数f(x,y)。%a,b分别为取值范围的左右端点。 %ya为给定初值y(a)。%N为迭代步数。%h为步长。 %输出值为对应每个节点的近似值。 h=(b-a)/N; T=zeros(1,N+1); Y=zeros(1,N+1); T=a:h:b; Y(1)=ya; for j=1:N Y(j+1)=Y(j)+h*feval('f',T(j),Y(j)); End T=[T' Y']
称阶方法。当 p＝1 时，一阶泰勒公式就是欧拉公式。 p 阶方法局部截断误差为：
Tn 1 y ( xn ) hf ( xn , y ( xn )) h
p
h
2
f
(1)
2!
( xn , y( xn ))
f
( p 1)
p!
( xn , y ( xn )) y ( xn 1 )
xn
欧拉方法 y 二阶泰勒 y 三阶泰勒 y 四阶泰勒 y
n n n
n
精确解 1 0.003322 1.013159 1.029143 1.050718 1.077217 1.107932 1.142165 1.179274 1.218689 1.259921
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
第六章 常微分方程初值问题初步
6.1 基本理论与Euler方法 6.2 Runge－Kutta方法（龙格－库塔法） 6.3用Matlab求解常微分方程的经典例 6.4常微分方程的解析解 6.5 方程组与高阶问题
6.1 基本理论与Euler方法
6.1.1 常微分方程的初值问题 考虑一阶微分方程的初值问题
在Euler算法程序的最后改为 Y1=sqrt(1+2*T); T=[T' Y' Y1']; 可得到解析解和迭代法计算结果比较
6.1.3 微分方程的几何解释 在 xy 平面上，积分曲线上任一点(x,y)的切线斜率等于函数 f(x,y)的值， 如果按函数 f(x,y)在 xy 平面上建立一个向量场，则积分曲线上每一点的切线 方向与向量场在该点的方向相一致。 从初始点 P0 ( x0 , y0 ) 出发，先依向量场在该点的方向推进到 x x1 上一点
泰勒级数法的 MATLAB 实现（四阶泰勒法）程序 通过计算 y 、 y 和 y ，并在每一步中使用泰勒级数，求区间[a , b]内的 y f ( x , y ) 初值问题： 。 y ( a ) y0 function T4=taylor4(df,a,b,ya,N)
(3) (4)
在牛顿-柯特斯公式中取 n 1 得梯形公式为：
2x y' y (0 x 1) y 【例 6.1】 利用程序解决初中问题： 。 y (0) 1 建立一个 f.m 的文件，
function z=f(x,y) z=y-2*x/y; 取迭代步骤为 10，在命令窗输入 euler('f',0,1,1,10) 实际上，该问题有解 y 似计算和精确值进行比较： 1 2 x （通过变量分离法可求出） ，对近
n
偏离了原来的积分曲线，可见欧拉方法的确精度不高。 局部截断误差： （利用泰勒展开式可得）
Tn 1 y ( xn 1 ) yn 1
( xn xn 1 )
h
2
y( )
h
2
2
2
y( xn )
假定 y(x)区间[a,b]上允分光滑，并令 M max y( xn )
dy f (t , y ) dt y( t ) y 0 0
如果存在实数 L 0 ，使得 | f ( x , y1 ) f ( x , y2 ) | L | y1 y2 | 则称 f 关于 y 满足 Lipschitz 条件，其中 L 称为 Lipschitz 常数。
(3)
4 x ( y ) y
3
2
y (4) y(3)
2 y
2
( xy
3 y )
12 y y
3
( xy y )
12 x ( y ) y
4
3
取 p 4 ，利用四阶泰勒公式，取 0.1，计算结果如下 表, 用四阶泰勒公式计算得结果 xn yn y ( xn ) 0.1 0.2 0.3 1.0954 1.1832 1.2649 1.0954 1.1832 1.2649
与欧拉方法相比，四阶泰勒公式精度较高。
取 p 4 ，利用四阶泰勒公式，取 0.1，计算结果如下 【例】 用泰勒方法解初值问题： 2 x y 2 3 y y (0) 1 解 可以求出精确解 y ( x )
3
1 xn 。
2
用欧拉方法、二阶、三阶、四阶泰勒方法计算的结果及精确值
6.1.7 梯形方法的求解 梯形方法是隐式的，用迭代求解，用欧拉方法提供的迭代初值的迭 代公式为：
y (0) y hf ( x , y ) n n n n1 ( k 1) h (k ) yn 1 yn [ f ( x n , yn ) f ( x n 1 , y n 1 ] 2 ( k 0,1, 2,)
定 理 6.1 设 f 在 区 域 D {( x , y ) | a x b, y R } 上 连 续 ， 关 于 y 满 足 Lipschitz 条件，则对任意的 x0 [a , b] ， y0 R ,常微分方程的初值问题在
x [a , b]存在连续可微解。
6.1.2 Euler 基本思想 可用差商代替，所以 dt y ( t 0 h) y ( t 0 ) f ( t 0 , y0 ) h 若令 y1 y0 hf ( t 0 , y0 ) ，而且可以得到迭代格式 yn 1 yn hf ( t n , yn ) 其中的 yn 就是 y ( t n ) 。 由数值微分的基本思想，
1 1.003333 1.013234 1.029346 1.051119 1.077866 1.108886 1.143468 1.180960 1.220785 1.262448
1 0.003322 1.013202 1.029241 1.050896 1.077498 1.108334 1.142704 1.179962 1.219683 1.260938
6.2.2 龙格-库塔方法的基本思想 利用 f(x,y) 在某些点处的值的线性组合构造公式， 使其按泰勒展开后与初 值问题的解的泰勒展开比较， 有尽可能多的项完全相同以确定其中的参数， 从而保证算式有较高的精度。 平均斜率 考察均差
此迭代格式是收敛的。
6．2
龙格-库塔方法
h
2 (1) fn
6.2.1 泰勒级数法 yn 1 yn hf n h
p
2!
p!
fn
( p 1)
，其中 f n f ( xn , yn ) ，
fn
(k )
dy f ( x, y) (k ) f ( xn , yn ) 称为初值问题 dx 的数值解的泰勒公式，又 y( x ) y 0 0
1 1.00867 1.019824 1.039054 1.063754 1.093211 1.123781 1.163444 1.202845 1.244311 1.287372
1 1.003333 1.013180 1.029171 1.050751 1.077252 1.107965 1.142194 1.179297 1.218706 1.259930
2 梯形方法的几何意义
用几何直观来说明梯形方法的平均化思想。
6.1,6 梯形方法的精度 上图中，设顶点 Pn ( xn , yn ) 在积分曲线 y=y(x)上，用欧拉公式过点 Pn 以斜率 f ( xn , yn ) 引折线交 x=xn+1 得顶点 A； 而后退的欧拉方法则以 点 Q(xn+1,yn+1)的斜率 f ( xn , yn ) 从顶点 Pn 引折线交 x 得另一顶点 B， 从图中可以看出， 和 B 两点均偏离点 Q 比较远， A 然而 AB 的中点 Pn+1 却相当接近点 Q，梯形方法改善了精度。
a xb
则 Tn 1 M
h
2
o( h ) ， 即欧拉方法的局部截断误差与 h2 同阶。
2
2
6.1.5 后退的欧拉公式
yn 1 yn hf ( xn 1 , yn 1 ) 后退欧拉公式： y ( x 0 ) y0
其局部截断误差：Tn 1 o( h ) 后退欧拉公式的求解
提高泰勒公式的阶 p，就可以提高计算结果的精度。
用泰勒公式求解初值问题：
2x y' y (0 x 1) y y (0) 1
解 直接求导得 2x y y y 2 y y 2 ( x xy ) y
y
(3)
y
ቤተ መጻሕፍቲ ባይዱ
2 y
2
( xy 2 y )
2
( n 0,1,)
欧拉公式是关于 yn 1 的一个直接的计算公式，这种公式称为显式 的；而后退欧拉公式是隐式的。 根本方法：用迭代法将隐式方程逐步显式化，迭代格式为：
yn 1
( k 1)
yn hf ( xn 1 , yn 1 ),
(k )
( k 0,1,)
后退欧拉公式的精度，后退欧拉公式的局部截断误差： 假设 yn y ( xn )
2
y ( x n 1 ) yn 1
h
2
y( xn )
6.1.5 梯形公式 梯形公式： yn 1 yn
2
h
2 实质：将欧拉公式与后退欧拉公式的误差估计式算术平均，消除误
差的主要部分
h yn
[ f ( xn , yn ) f ( xn 1 , yn 1 )]后退欧
P1 ， 再从 P1 依向量场方向推进到 x x2 上一点 P2 ， 这样下去便可作出一条折线
P0 P1 P2 。设已作出折线的顶点为 Pn ，过 Pn ( xn , yn ) 依向量场的方向再推进到
Pn 1 ( xn 1 , yn 1 ) 。
yn 1 yn x n 1 yn
欧 步长
f ( x n , yn )
拉 公

式 ：
yn 1 yn hf ( xn , yn ) ， 其中 h xi xi 1为
可 逐 步 y1 y0 hf ( x0 , y0 )
y2 y1 hf ( x1 , y1 )
计
算
出
：



6.1.4 欧拉公式的精度 假设 yn y ( xn ) ，即顶点 Pn 落在积分曲线 y y ( x ) 上，按欧拉方 法作出的折线 Pn Pn 1 便是 y y ( x ) 过点 P 的切线，定出的顶点 Pn 1 明显
%df 为 y 的一阶到四阶微商序列=[y' y'' y''' y'''']。
function T4=taylor4(df,a,b,ya,N) %df为y的一阶到四阶微商序列=[y' y'' y''' y'''']。%y'=f(x,y)。%a,b左右端点。%N 为迭代步数。%h为步长。%ya为初值。 h=(b-a)/N; T=zeros(1,N+1); Y=zeros(1,N+1); T=a:h:b; Y(1)=ya; for j=1:N D=feval(df,T(j),Y(j)); Y(j+1)=Y(j)+h*(D(1)+h*(D(2)/2+h*(D(3)/6+h*D(4)/24))); end T4=[T' Y']
【例】 利用四阶泰勒法求解初值问题
1 y ( x y ) 2 y (0) 1
其中[a , b] [0,1]。 解 建立 M 文件 df.m function z=df(x,y) z=[(x-y)/2 (2-x+y)/4 (-2+x-y)/8 (2-x+y)/16]; 在 MATLAB 运行窗口输入 taylor4('df',0,1,1,10)