微分方程数值解法(戴嘉尊_第二版)习题讲解
成都信息工程学院>>精品课程>>微分方程数值解
微分方程数值解
习题解答
杨韧 吴世良(编)
成都信息工程学院
数学学院
二O一O年四月编写
电子文档制作：成都信息工程学院 数学学院 杨韧 吴世良，2010 年 4 月
成都信息工程学院>>精品课程>>微分方程数值解
目 录
第一章 常微分方程数值解 ......................................................................3
第二章 抛物型方程的差分方法 ..............................................................8
第三章 椭圆型方程的差分方法 ............................................................16
第四章 双曲型方程的差分方法 ............................................................25
电子文档制作：成都信息工程学院 数学学院 杨韧 吴世良，2010 年 4 月
成都信息工程学院>>精品课程>>微分方程数值解
第一章 常微分方程数值解
1．解: 由欧拉公式得
yn1 yn hf (xn, yn) yn h( 2 2yn 2 ) yn 0.2y
n
2
1
1x
0.1
1x
n2
n
由梯形公式得
yn1 yn 2 yn 2 h [( yn hyn h[ f (xn, yn)f (xn1, yn1)]
1
2
2yn 2 )(1x 2y
2
n1
2
n1
)]
1 1 1x 1
n
2 hy 2
n1
1 2 h( 1 1x 2 1x 1
)
2
n
n1
hy 2 n1 yn1 yn hyn 2 1 2 h( 1 1x 2 1x 1 )
2
n
n1
114h(yn hy
n
2
h( 2 ))
1 2 1 1x 1
2
y
n1
1 x
n1
n
2h
欧拉公式计算结果
x
n
y
n
y(xn )
y(xn)y
n
0 0 0 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1000 0.1970 0.2854 0.3609 0.4210 0.4656 0.4957 0.5137 0.5219 0.5227 0.0990 0.1923 0.2752 0.3448 0.4000 0.4412 0.4698 0.4878 0.4972 0.5000 0.0010
0.0047
0.0102
0.0160
0.0210
0.0244
0.0259
0.0259
0.0247
0.0227
梯形公式计算结果
x
n
y
n
y(xn )
y(xn)y
n
0 0 0 0
电子文档制作：成都信息工程学院 数学学院 杨韧 吴世良，2010 年 4 月
微分方程数值解法
微分方程数值解法微分方程数值解法微分方程数值解法【1】摘要：本文结合数例详细阐述了最基本的解决常微分方程初值问题的数值法，即Euler方法、改进Euler法，并进行了对比，总结了它们各自的优点和缺点，为我们深入探究微分方程的其他解法打下了坚实的基础。
关键词:常微分方程数值解法 Euler方法改进Euler法1、Euler方法由微分方程的相关概念可知，初值问题的解就是一条过点的积分曲线，并且在该曲线上任一点处的切线斜率等于函数的值。
根据数值解法的基本思想，我们取等距节点，其中h为步长，在点处，以为斜率作直线交直线于点。
如果步长比较小，那么所作直线与曲线的偏差不会太大，所以可用的近似值，即：，再从点出发，以为斜率作直线，作为的近似值，即：重复上述步骤，就能逐步求出准确解在各节点处的近似值。
一般地，若为的近似值，则过点以为斜率的直线为：从而的近似值为：此公式就是Euler公式。
因为Euler方法的思想是用折线近似代替曲线，所以Euler方法又称Euler折线法。
Euler方法是初值问题数值解中最简单的一种方法，由于它的精度不高，当步数增多时，由于误差的积累，用Euler方法作出的折线可能会越来越偏离曲线。
举例说明：解： ,精确解为：1.2 -0.96 -1 0.041.4 -0.84 -0.933 0.9331.6 -0.64 -0.8 0.161.8 -0.36 -0.6 0.242.0 0 -0.333 0.332.2 0.44 0 0.44通过上表可以比较明显地看出误差随着计算在积累。
2、改进Euler法方法构造在常微分方程初值问题 ,对其从到进行定积分得:用梯形公式将右端的定积分进行近似计算得：用和来分别代替和得计算格式:这就是改进的Euler法。
解：解得：由于 ,是线形函数可以从隐式格式中解出问题的精确解是误差0.2 2.421403 2.422222 0.000813 0.021400.4 2.891825 2.893827 0.00200 0.051830.6 3.422119 3.425789 0.00367 0.094112.0 10.38906 10.43878 0.04872 1.1973通过比较上表的第四列与第五列就能非常明显看出改进Euler方法精度比Euler方法精度高。
微分方程数值解法答案
微分⽅程数值解法答案包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。
解答问题关键在过程，能够显⽰出你已经掌握了书上的内容，知道了解题⽅法。
这次考试题⽬的类型：20分的选择题，主要是基本概念的理解，后⾯有五个⼤题，包括差分格式的构造、截断误差和稳定性。
习题⼀1．略2． y y x f -=),(，梯形公式：n n n n n n y hh y y y h y y )121(),(2111+-+=+-=+++，所以0122)1(01])121[()121()121(y hh y h h y h h y hhn h h n n n +--+--+-+=+-+==+-+= ，当0→h 时，x n e y -→。
同理可以证明预报-校正法收敛到微分⽅程的解.3．局部截断误差的推导同欧拉公式;整体截断误差:++++++-++≤1),())(,(11111n nx x n n n n n n n dx y x f x y x f R εε11)(++-++≤n n n y x y Lh R ε，这⾥R R n ≤ ⽽111)(+++-=n n n y x y ε，所以 R Lh n n +=-+εε1)1(，不妨设1()]11111[1111101---++-+-+-≤≤-+-=n n n n Lh Lh Lh R Lh Lh R Lh εεε ]1[2)(02)(00-+≤--x X L x X L eLh R eε4．中点公式的局部截断误差： dx x y x f hx y h x f x y x f yx y n n x x n n n n n n))](,(2)(,2())(,([)(11*1?+++-=-++dx x y x f hx y h x f h x y h x f h x y x y dxx y x f hx y h x f hx y h x f h x y h x f x y x f n n n n x x n n n n n n n x x n n n n n n n n))](,(2)(,2())2(,2([)]2()([))](,(2)(,2())2(,2())2(,2())(,([11++-++++'-'=++-+++++-=??++所以上式为+--+''=?++dx hx x x y e n nx x n n n )2()(11θdx x y x f h x y h x f h x y h x f n n n n x x n n n n))](,(2)(,2())2(,2([1++-++?+ 3218)(LMh h x y Lh e n n ≤+''≤+?中点公式的整体截断误差：dx y x f hy h x f x y x f y x y y x y n n x x n n n n n n n n)],(2,2())(,([)()(111?+++-+-=-++dxy x f hy h x f x y x f h x y h x f x y x f hx y h x f x y x f y x y n n n n n n n n x x n n n n n n n n))],(2,2()))(,(2)(,2()))(,(2)(,2())(,([)(1++-+++++-+-=?+因⽽n n n L h Lh R εεε)21(1+++≤+，R L h Lh n n +++≤-122)21(εε≤≤])21()21(1[2)21(1222222022-+++++++--+++n nL h Lh L h Lh Lh Lh RL h Lh ε )1(00-+≤--x X L x X L e LhR eε 5．略 6．略 7．略8．（1）欧拉法：2.0≤h ；四阶Runge-Kutta ⽅法：278.0≤h （2）欧拉法：3 54≤h ；四阶Runge-Kutta ⽅法：3556.5≤h（3）欧拉法：1≤h ；四阶Runge-Kutta ⽅法：278.0≤h 9．略 10．略习题21．略 2．略 3．略4．差分格式写成矩阵形式为：n n M n M n n n M n M n n e u u u u r t r r r t r r r t r r r t u u u u +?--------= --+-+-++12211221121212121 αβαααβαααβαααβ矩阵的特征值为：)cos(221Mj r r t j πααβλ+-?-=，要使格式稳定，则特征值须满⾜t c j ?+≤1λ，即21≤r α5．利⽤泰勒展式可以得到古典隐式差分格式的截断误差为)(2h t O +?。
微积分(下册)第二版 第9章 微分方程
又y 0是原方程的解 ,
方程的解为: y C 1 x2 , C为任意常数.
二、齐次方程 形式：
dy dx
f
y x
解法 作变量代换 u y , 即 y xu, x
dy u x du , u x du f (u),
dx
dx
dx
即
1 du f (u) u
1 dx. x
du f (u) u
◆微分方程的分类： 分类1: 常微分方程, 偏微分方程.
分类2: 一阶微分方程 F ( x, y, y) 0, y f ( x, y); 高阶微分方程 F ( x, y, y,, y(n) ) 0,
y(n) f ( x, y, y,, y(n1) ).
2023/8/29
2
分类3: 线性与非线性微分方程: y P( x) y Q( x), x( y)2 2 yy x y 0.
1 dx, x
例3
求解 dy x y 1. , 令u y，
dx
x
x
则y xu， dy u x du ,
dx
dx
u x du 2 u, dx
即 x du 2, dx
x du 2, dx
du 2 dx, 两边积分,得 u ln x2 c, x
如: y y, 通解 y cex;
y y 0, 通解 y c1 sin x c2 cos x.
(2)特解: 通解中任意常数为确定值的解.
初始条件: 用来确定任意常数的条件.
2023/8/29
5
§9.2 可分离变量的微分方程 一、 可分离变量的微分方程 二、 齐次方程
一、可分离变量的微分方程
dx
dx
2 y3 dy 2xy2 2x3 dx
微分方程数值解
第四章 微分方程数值解4.1 算 法当常微分方程能解析求解时，可利用Matlab 符号工具箱中的功能找到精确解. 见下例求解方程20y y y '''++=. 键入：syms x y %定义符号变量 diff_equ= ‘D2y+2*Dy -y=0’; %D2y 表示,y ''Dy=y 'y=dsolve (diff_equ, ‘x’)%定义x 为自变量y=cl*exp ((2^(1/2)-1)*x+c2*exp (-(2^(1/2)+1)*x) %表达式中含c1与c2，表示通解.%初始条件为y (0)=0,y '(0)=1时，按如下方式调用y=dsolve (diff_equ, ‘y (0)=0’, ‘Dy (0)=1’, ‘x’)y=1/4*2^(1/2)*exp ((2^(1/2)-1)*x)— 1/4*2^(1/2)*exp (-(2^(1/2)+1)*x) %画出函数y=y (x)的图形ezplot (y,[-2,2])图形具体形式请上机试之.在方程无法获得解析解的情况下，可方便地获得数值解. 下面的例子说明用Matlab 求数值解的方法及应注意的问题.[例1] 求解范德堡（vander pol ）方程222(1)0d x dx x x dt dtμ--+=求解高阶方程，必须等价地变换为一阶微分方程组，对本例，通过定义两个新的变量，实现这一变换1,2/y x y dx dt ==则令 1/2dy dt y =22/(11)*21dy dt y y y μ=--编写求解程序分为两部分，第一部分为待求解的方程，存盘的文件名为<待求解方程的函数名.m >，第二部分为求解主程序，本例中取名为main1.m.首先编写待求解方程的M -文件. 文件存盘名为“vdpol.m ”. function yprime=vdpol (,)t y yprime (1)=y (2);mu=2yprime (2)=mu*(1(1)^2)*(2)(1)y y y --;yprime=[yprime (1);yprime (2)];说明 函数yprime=vdpol (,)t y 中. 定义t 为自变量，y 的形式取决于求解方程的阶数，本例中，[(1),(2)],(1)y y y y =为解向量，(2)y 为导数向量. yprime (1)(1)y '=， yprime (2)(1)y ''=，函数返回vander pol 方程的导数列向量. 因为所求结果为方程数值解，所以各向量维数只有在主程序求解时定下精度后才能确定.主程序定名为main1.m ，你可用你所喜欢的其它名子，但vdpol.m 除外. clear functions%调试程序时，放置这一语句是必要的. 它清除前边已编译的存在于内存中的废弃程序 [,t y ]=ode23 (‘vdpol’,[0,30],[1,0]); y1=y (:,1); %解曲线. y2=y (:,2); %解曲线的导数. polt (,1,,2,t y t y ‘_ _’)说明 龙格_库塔的2阶与4阶改进型求解公式的实现，其指令分别为： [,t x ]=ode23 (‘f ’,,0,ts x options) [,t x ]=ode45 (‘f ’,,0,ts x options)其中ts 可由系统依据精度要求自动设定，亦可由使用者依据实际需要自己确定，分别说明之. （1）若令[0,1,,]ts t t tf =，则输出在指定时刻0,1,,t t tf 给出，当0::ts t k tf =时，输出在区间[0,]t tf 的等分点上给出，k 为步长.（2）若[0,],0ts t tf t =为自变量初值，tf 为终值，此时，options 决定自变量t 的维数，t 中的时间点不是等间隔的，这是为了保证所需的相对精度，积分算法改变了步长. 用于设定误差限的参数options 可缺省，此时系统设定相对误差为310-，绝对误差为610-，若自行设定误差限，可用如下语句： options=odeset (‘reltol’,rt , ‘abstol’,at ) 这里的rt 与at 分别为设定的相对与绝对误差.须注意的是无论用哪种方法确定t 的取值方式，0,t tf 必须由使用者确定且应与0x 相匹配. 0x 为初始条件，本例中0[1,0]y =，因为[0,30]ts =，这意味着解曲线(0)1y =，(0)0y '=，一般说，当解n 个未知函数的方程组时，0x 为n 维向量，共含有n 个初始条件.两个输出参数是列向量t 与矩阵x ，它们具有相同的行数，而矩阵x 的列数等于方程组的个数，本例中y 的列数为2，其中，(:,1)y 为自变量t 上各点函数值，(:,2)y 为t 上各点导数值.最后，提请读者注意的是：ode45也不总是比ode23好，在很多时候，低阶算法更有效，有关微分方程数值解法的更进一步信息，请参考数值分析方面的书籍. 有些参考书提供了一些关于算法选择和如何处理那些时间常数变化范围大的病态方程的非常实用的算法.4.2 欧拉与龙格-库塔方法设有一阶方程与初始条件0(,)()y f x y y x y '=⎧⎨=⎩ （4.1） 其中(,)f x y 适当光滑，关于y 满足Lipschitz 条件，即存在L 使1212(,)(,)f x y f x y L y y -≤-则（4.1）式的解存在且惟一.关于()y y x =的解析解一般难以求到或根本无解析解，因而，实际问题中，通常，采用差分的方法. 在一系列离散点12n x x x <<<<上寻求其数值近似解12,,,,n y y y .相邻两个节点间的间距1n n n h x x +=-称为步长，一般地取等步长h ，则0,1,2,n x x nh n =+=1、欧拉方法在区间1[,]n n x x +上用差商1()()n n y x y x h+-代替（4.1）式中y '，对(,)f x y 中x 在1[,]n n x x +上取值n x 还是1n x +，而形成向前欧拉公式与向后欧拉公式.（1）向前欧拉公式(,)f x y 取左端点n x ，得如下公式1()()(,())n n n n y x y x hf x y x +≈+ （4.2）从0x 点出发，由初值00()y x y =代入（4.2）求得1000(,)y y hf x y =+ （4.3）反复利用（4.2），有1(,) 0,1,2,n n n n y y hf x y n +=+= （4.4）其几何意义如图4.1所示.图中()y y x =为方程（4.1）的精确 解曲线，其上任意点(,)x y 处切线斜率为(,)f x y . 从初值000(,)P x y 点出发，用该点斜率00(,)f x y 作一直线段，在1x x = 处得到111(,)P x y ，1y 由（4.2）式确定， 再从111(,)P x y 出发，以11(,)f x y 为斜率 作直线段，在2x x =处得到222(,)P x y ，y0yO0x 误差4x1x 2x 3x x()y y x =3()y x32()y y x - 0P1P2P3P 4P3y这一过程继续下去，形成折线012P PP ，作为积分曲线()y y x =的近似，用1()n y x + 图4.1 表示在1n x +处的精确值，1n y +为解的近似值，不难得到23211()()()()2n n n h y x y y x O h O h ++''-=+=这一误差称为局部截断误差. 若一种算法局部截断误差为1()P O h +，则称该算法具有P 阶精度，所以向前欧拉公式具有1阶精度.（2）向后欧拉公式若(,)f x y 中x 取1[,]n n x x +中的1n x +，则有如下公式：111(,) 0,1,2,n n n n y y hf x y n +++=+= （4.5）称式（4.5）为向后欧拉公式，因为此式中1n y +未知，故称其为隐式公式，无法用其直接计算1n y +，一般用向前欧拉公式产生初值.(0)11(,) 0,1,2,n n n n y y hf x y n ++=+=再按下式迭代(1)()111(,),0,1,,0,1,k k n n n n y y hf x y k n ++++=+==其误差估计如下23211()()()()2n n n h y x y y x o h o h ++''-=+=精度亦为1阶，将向前欧拉公式（4.4）与向后欧拉公式（4.5）及它们的误差的几何说明作一对比，是十分有益的，见图4.2. 为讨论局部截断误差，在图4.2中设点 (,)n n n P x y 落在积分曲线()y y x =上，按式 （4.4）及式（4.5）分别得1n P +点为A 与B ， 且,A B 点一定在积分曲线()y y x =上相应 点θ的上、下两边,所以将式（4.4）与（4.5） 平均之，一定能得到更好的结果. （3）梯形公式 将向前与向后欧拉公式加以平均得到所 图4.2谓梯形公式111[(,)(,)] 0,1,2,2n n n n n n hy y f x y f x y n +++=++= （4.6）其局部截断误差为3()O h ，具有2阶精度.（4）改进的欧拉公式为使计算简单，又免去迭代的繁复，将公式（4.6）简化为两步1(,)n n n n y y hf x y +=+x 1n x +n x y()y y x =n PAθB111[(,)(,)], 0,1,2,2n n n n n n hy y f x y f x y n +++=++= （4.7）或写为1121211()2(,)(,)n n n n n n h y y k k k f x y k f x y hk ++⎧=++⎪⎪=⎨⎪=+⎪⎩0,1,2,n= （4.8）最后指出，上述欧拉方法可推广至微分方程组，如0000(,,)(,,)()()y f x y z z g x y z y x y z x z '=⎧⎪'=⎪⎨=⎪⎪=⎩ 向前欧拉公式为11(,,)(,,)n n n n n n n n n n y y hf x y z z z hg x y z ++=+⎧⎨=+⎩ 0,1,2,n=2、龙格_库塔方法 由微分中值定理1[()()]/(),01n n n y x y x h y x h θθ+'-=+<<又因为(,)y f x y '=，所以()(,())n n n y x h f x h y x h θθθ'+=++从而有1()()(),()n n n n y x y x hf x h y x h θθ+=+++ （4.9）令(,())n n k f x h y x h θθ=++，称其为区间1[,]n n x x +上的平均斜率，由（4.9）可知，给出一种平均斜率，可相应导出一种算法. 向前欧拉公式中(,)n n k f x y =，精度低. 改进欧拉公式中取111((,)(,))2n n n n k f x y f x y ++=+，精度提高，下面，我们在区间1[,]n n x x +内多取几个点，将其斜率加权平均，就能构造出精度更高的计算公式，公式的推导不再具体给出，只开列具体结果. （1）2阶龙格_库塔公式11122121()(,)(,),0,1n n n n n n y y h k k k f x y k f x ah y hk a λλββ+=++⎧⎪=⎨⎪=++<<⎩ （4.10） 其中12211,,12a aβλλλ+===，由于4个未知数只有3个方程，所以解不惟一，若令121,12a λλβ+===，即得改进的欧拉公式，具有2阶精度.（3）4阶龙格_库塔公式 只给出精典格式中最常用的一种.112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y k k k k k f x y h h k f x y k h h k f x y k k f x h y hk +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩ （4.11） 其计算精度为4阶4.3 模型与实验1、模型与问题 [例2] 单摆运动图4.3中一根长l 的细线，一端固定，另一端悬挂质量为m 的小球，在重力作用下，小球处于竖直的平衡位置. 现使小球偏离平衡位置一个小的角度θ，然后使其自由运动，在不考虑空气阻力情形下，小球将沿弧线作周期一定的简谐运动. 0θ=为平衡位置，在小球摆动过程中，当与平衡位置夹 角为θ时，小球所受重力在其动运轨迹的分量为sin mg θ- （负号表示力的方向使θ减少），由牛顿第二定律可得微分方 程()sin ml t mg θθ''=- （4.12） 设小球初始偏离角度为0θ，且初速为0，式（4.12）的初 始条件为0(0),(0)0θθθ'== （4.13）当0θ不大时，sin θθ≈，式（4.12）化为线性常系数微分方程 图4.30g lθθ''+= （4.14）解得θlθmg0()t θθ= （4.15）简谐运动的周期为2T =. 现在的问题是：当0θ较大时，仍用θ近似sin θ，误差太大，式（4.12）又无解析解，试用数值方法在0030,10θθ==两种情况下求解，画出()t θ的图形，与近似解（4.15）比较，这里设25cm l =. [例3] 捕食与被捕食当鲨鱼捕食小鱼，简记为乙捕食甲，在时刻t ，小鱼的数量为()x t ，鲨鱼的数量为()y t ，当甲独立生存时它的（相对）增长率与种群数量成正比，即有()()x t rx t '=，r 为增长率，而乙的存在使甲的增长率r 减少，设减少率与乙的数量成正比，而得微分方程()()(())x t x t r ay t rx axy '=-=- （4.16）比例系数a 反映捕食者掠取食饵的能力.乙离开甲无法生存，设乙独自存在时死亡率为d ，()()y t dy t '=-，甲为乙提供食物，使乙的死亡率d 降低，而促其数量增长，这一作用与甲的数量成正比，于是()y t 满足()(())y t y d bx t dy bxy '=-+=-+ （4.17）比例系数b 反映甲对乙的供养能力，设若甲，乙的初始数量分别为00(0);(0)x x y y == （4.18） 则微分方程（4.16），（4.17）及初始条件（4.18）确定了甲，乙数量()x t 、()y t 随时间变化而演变的过程，但该方程无解析解，试用数值解讨论以下问题：（1）设001,0.5,0.1,0.02,25,2r d a b x y ======，求方程（4.16），（4.17）在条件（4.18）下的数值解，画出(),()x t y t 的图形及相图(,)x y ，观察解(),()x t y t 的周期变化，近似确定解的周期和,x y 的最大、小值，近似计算,x y 在一个周期内的平均值. （2）从式（4.16）和（4.17）消去dt 得到()()dx x r ay dy y d bx -=-+ （4.19） 解方程（4.19），得到的解即为相轨线，说明这是封闭曲线，即解确为周期函数. （3）将方程（4.17）改写为1()[]y x t d b y'=+ （4.20）在一个周期内积分，得到()x t 一周期内的平均值，类似可得()y t 一周期内的平均值，将近似计算的结果与理论值比较. 2、实验（1）方程（单摆问题）0sin (0),(0)0ml mg θθθθθ''=-⎧⎨'==⎩无解析解，为求其数值解，先将其化为等价的一阶方程组. 令12,x x θθ'==，方程化为1221102sin (0),(0)0x x g x x l x x θ'=⎧⎪⎪'=-⎨⎪==⎪⎩ 其中，09.8,25,g l θ==为100.1745≈（弧度）与300.5236≈（弧度）两种情况，具体编程如下：先以danbai.m 为文件名存放待解方程. 键入： function xdot =danbai (t,x) %x=[x (1),x (2)],g=9.8;1=25; %x (1)为解向量, x (2)是导数 xdot (1)=x (2); %xdot (1)=x '(1)=θ' xdot (2)=-g/1*sin (x (1)); %xdot (2)=θ''xdot=[xdot (1);xdot (2)]; %必须将导数向量以列向量形式给出.再以主程序(求数值解)调用待求方程，主程序用main2.m 为文件名存盘，其代码如下： clear functions[t,x]=ode23 (‘danbai’,[0,10],[0.1745,0])%只计算0100.1745θ==（弧度）的情形. %对近似解0()cos ,t wt w θθ==2T = w=sqrt (9.8/25); y=0.1745*cos (w*t);[t,x (:,1),y] %显示数据，无分号. hold on %欲在同一幅图中画近似解. plot (t,x (:,1), ‘b’) %画数值解，绿色 plot (t,y, ‘g*’) %用*号，红色画近似解. Hold off（2）食饵_捕食者 方程()x t rx axy '=- ()y t dy bxy '=-+可化为如下形式00x r ay x d bx y y ⎡⎤-⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥-+⎣⎦⎣⎦⎣⎦初始条件00(0),(0)x x y y ==表示为00(0)(0)x x y y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦以shier.m 存盘如下代码 function xdot=shier (t,x) r=1;d=0.5;a=0.1;b=0.02;xdot=diag ([r-a*x (2),-d+b*x (1)])*x;%x=(1)(2)x x x y ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦,xdot=x y ⎡⎤⎢⎥⎢⎥⎣⎦以main3.m 存盘如下代码. clear functions ts=0:0.1:15; x0=[25,2];[t,x]=ode45 (‘shier’,ts,x0);[t,x] %显示数据t,x,y plot (t,x)gird %加网格线 gtext (‘x(t)’),gtext (‘y(t)’), %用点鼠标方式pause,figure (2) %将x 1(t),x 2(t)放至指定点 plot (x (:,1),x (:,2)) %以x (1)与x (2)为坐标点画相图xlabel (‘x’),ylabel (‘y’)可以猜测，1()x t 与2()x t 是周期函数，相图12(,)x x 是封闭曲线，从数值解可近似定出周期为110.7,x 的最大和最小值分别为99.3与22.0,x 的最大和最小值分别为28.4和2.0，为求1x 与2x 在一个周期的平均值，只需键入： y1=x (1:108,1); %x 1周期为10.7. x1p=trapz (y1)*0.1/10.7, %trapz (y1)返回按 y2=x (1:108,2); %梯形法对y1的积 x2p=trapz (y2)*0.1/10.7, %分值10711()1(1)2i i y i y i x =++∆∑ 可得1225,10x x == 对方程()()dx x r ay dy y d bx -=-+化为d bx r aydx dy x y-+-=积分得解()()d bx r ay x e y e c --= （4.21）即为原方程组的相轨线，其中c 由初始条件确定. 为说明上述相轨线是封闭的，令();()d bx r ay f x x e g y y e --==设其最大值分别为,m m f g ，若00,x y 满足00(),()m m f x f f y g ==则有00,d rx y b a==（令0,0f g ''==，可解出00,x y ），又当m m f g c ⋅≥时，相轨线（4.21）有意义. 当m m f g c =时，相轨线退化为一个点00(,)P x y ，对0m m c f g <<时，相轨线如图 4.4（c ），而图（a ），（b ）分别为()f x 与()g y 的图像，下面讨论如何由（a ），（b ）画（c ）.（a ） （b ） （c ）图4.4设(0)m m c pg p f =<<，若令0y y =，则有()f x p =，由（a ）知，12,x x ∃，使12()()f x f x p ==，且102x x x <<于是相轨线应过110(,)q x y 与220(,)q x y ，对12(,)x x x ∈，()f x p >，由()()()m m f x g x pg g y g =⇒<，令()g y q =，又由（b ）知，存在12,y y 使12()()g y g y q ==，且102y y y <<，于是相轨线又过31(,)q x y 与42(,)q x y 两点，所以对12,q q 间每个x ，轨线总要通过纵坐标为12102,()y y y y y <<的两点，于是相轨线是一条封闭曲线.相轨线封闭等价于(),()x t y t 是周期函数，记周期为T ，为求其在一个周期内的平均值(,)x y ，由1()()yx t d b y=+两边在一个周期内积分有((0)())y y T =：011ln ()ln (0)()T y T y dT d x x t dt T T b b b-⎡⎤==+=⎢⎥⎣⎦⎰ 同理()f xmfP1x 0x 2xxyy 1y ()()n m n mf x fg y g ==2yyqm g()g vx1y 0y 2y1xx 2x1q 3q2q00(,)P x yr y a=从而 00,x x y y ==即(),()x t y t 的平均值恰为相轨线中心点p 的坐标，提请读者注意的是：这里的0x 与0y 与初始条件中的00,x y 不是一件事.注意到,,,r d a b 在生态学上的意义，上述结果表明，捕食者的数量与食饵的增长率成正比，与它掠取食饵的能力成反比，食饵的数量与捕食者死亡率成正比，与它供养捕食者的能力成反比.3、练习内容（1）编写改进欧拉公式求微分方程数值解的程序，并用其与ode23求下列微分方程数值解，对二者作出比较.a ）22,(0)0y x y y '=-=或(0)1y =.b)222()0x y xy x n y '''++-=2()2,()22y y πππ'==-（Bessel 方程，这里令12n =，其精确解为sin y =. c)cos 0,(0)1,(0)0y y x y y '''+===. （2）倒圆锥形容器，上底面直径为1.2m ，容器的高亦为1.2m ，在锥尖的地方开有一直径为3cm 的小孔，容器装满水后，下方小孔开启，由水利学知识可知当水面高度为h 时，水从小孔中流出的速度为v g =为重力加速度，若孔口收缩系数为0.6（即若一个面积单位的小孔向外出水时，水柱截面积为0.6），问水从小孔中流完需多少时间？2分钟时，水面高度是多少？（3）一只小船渡过宽为d 的河流，目标是起点A 正对着的另一岸上B 点，已知河水流速1v 与船在静水中的速度2v 之比为k .（a ）建立小船航线的方程，求其解析解.（b ）设12100m,1m/s,2m/s d v v ===，用数值解法求渡河所需时间，任意时刻小船的位置及航行曲线，作图并与解析解比较.（c ）若流速1v 为0,0.5,2 (m/s)，结果将如何？（4）研究种群竞争模型. 当甲、乙两个种群各自生存时，数量演变服从下面规律1212()(1),()(1)x y x t r x y t r y n n =-=- 其中，(),()x t y t 分别为t 时刻甲，乙两个种群的数量，12,r r 为其固有增长率，12,n n 为它们的最大容量，而当这两个种群在同一环境中生存时，由于乙消耗有限资源而对甲的增长产生影响，将甲的方程修改为1112(1)x y x r x s n n =-- （4.22） 这里1s 的含意是：对于供养甲的资源而言，单位数量乙（相对2n ）的消耗率为单位数量甲（相对1n ）消耗的1s 倍，类似地，甲的存在亦影响乙的增长，乙的方程应改为2212()(1)x y y t r y s n n =-- （4.23） 给定种群的初始值为 00(0),(0)x x y y == （4.24）及参数121212,,,,,r r s s n n 后，方程（4.22）与（4.23）确定了两种群的变化规律，因其解析解不存在，试用数值解法研究以下问题：（a ）设121212001,100,0.5,2,10r r n n s s x y ========，计算(),()x t y t ，画出它们的图形及相图(,)x y ，说明时间t 充分大以后(),()x t y t 的变化趋势（人们今天看到的已经是自然界长期演变的结果）.（b ）改变1200,,,r r x y ，但1s 与2s 不变（保持121,1s s <>），分析所得结果，若121.5(1),0.7(1)s s =>=<，再分析结果，由此你得到什么结论，请用各参数生态学上的含义作出解释.（c ）试验当120.8(1),0.7(1)s s =<=<时会有什么结果；当121.5(1), 1.7(1)s s =>=>时又会出现什么结果，能解释这些结果吗？。
第一章常微分方程初值问题数值解法
（1.2.3）
其中rn,k(t)为插值余项。 代到（1.2.2）式中得
u ( tn +1 ) = u ( tn ) +
舍去余项 并用uj代替u(tj)即得
∫
tn+1 tn
Ln , k ( t ) dt + ∫ t rn , k ( t ) dt
tn+1
n
(1.2.4 (1.2.5)
Rn , k = ∫
⎡ ∑ ⎣α u
j =0
j =0
αk ≠ 0
（1.2.1）
j n+ j
⎤ − hβ j f n + j ⎦ = 0（数值解满足的差分方程）
因此称（1.2.1）为多步法 或 k-步法。 又因为(1.2.1)关于 u n + j , f n + j 是线性的，所以称为线性多步法。 为使多步法的计算能够进行，除给定的初值u0 外，还要 知道附加初值u1,u2,…,uk-1 ，这可用其它方法计算。 若 β k = 0 则称(1.2.1)是显式的； 若 β k ≠ 0 则方法（1.2.1）是隐式的。 例如，一般线性二步法可写成：
f ( t , u ( t ) ) = Ln , k +1 ( t ) + rn , k +1 ( t )
其中rn,k+1(t)为插值余项。 同理即
un +1 = un + h ∑ bk +1i f ( tn −i +1 , un −i +1 )
i =0
k +1
其中
bk +1i
=∫ ∏
−1
j =0 j ≠i
0
k +1
20100319 第三讲：求导积分与微分方程数值解(2次课)
求导数运算的调用格式
[3] 参数方程求导 方法说明： 对参数方程x=x(t);y=y(t);先求出dy/dt和dx/dt 然后代入公式dy/dx= dy/dt / dx/dt 即可 应用示例： x t (1sin t ) 例6 求参数方程 y t cos t
syms t; x=t*(1-sin(t)); y=t*cos(t); ezplot(x,y); grid on; dx=diff(x,t); dy=diff(y,t); dydx=dy/dx; pretty(dydx) %下面在t=4.1处作出参数方程的切线（导数） hold on; t=4.1; x=eval(x); y=eval(y); plot(x,y,'ro'); k=eval(dydx); line([x,x+1],[y,y+k],'color','r')
18
程 数 值 解
第 三 讲 求 导 积 分 与 微 分 方
19
程 数 值 解
第 三 讲 求 导 积 分 与 微 分 方
That’s all~3Q!
求定积分运算的应用示例
例16 用两种方法求定积分
x=2:.1:5;y=log(x)./(x.^2); tt=trapz(x,y) 程 数 值 解 第 三 讲 求 导 积 分 与 微 分 方 %复合梯形公式
15
ln x / x 2 dx
2
5
fun=inline(' log(x)./(x.^2) ','x'); ss=quad(fun,2,5)
8Байду номын сангаас
z x3 y2 sin( xy),求3 z x3 例8 对
《微分方程数值解法》PPT课件
方程的解 U~n 。为了弄清差分格式（2.58）的稳定性条件， 给出稳定的定义：
对于任意给定的 0 ，存在与h, k 无关且依赖于 的
正数 ，使当
U~0 U 0 V 0
时，对于任何的 n0 nk T ，差分格式得到的解U~ n ,U n
满足不等式
U~n U n V n
连同初值条件：U
0 m
mk
, m
1,2,M
1
边值条件：U
n 0
U
n M
0, n
0,1,2,, N
逐层解出结点处的U 值。
现在对
h
ห้องสมุดไป่ตู้
，取二种 20
k
，使
r
k h2
5 和5 11 9
。图2.9
和图2.10中的曲线表示不同的时刻微分方程的精确解，图
中“ ”表示差分方程的解
（2.54）
下面我们先研究上式右边第二项，即差分方程的理论
解与计算机上解得的近似解之间的差别是随着n 的增大而
无限增大还是有所控制。如果这种差别是无限增加，则称
差分格式不稳定，显然不稳定的格式是不能使用的，因为
误差的无限增加淹没了真解。上例中的r 5 时就是差分
9
方程不稳定的情况。从差分方程，比如格式（2.29）可知，
如果差分方程为显式，则对所有的n ，An I ；如 An I
果 An
0
，
U U
n1 CnU
0
n
An1en , Cn
An1Bn
（2.58）
，则隐式格式可以写成显式形式。