整理微分方程详解.docx
《整理微分方程详解.docx》由会员分享，可在线阅读，更多相关《整理微分方程详解.docx（30页珍藏版）》请在冰豆网上搜索。
整理微分方程详解
第二章微分方程
本章学习目的：
本章的主要目的在于：
学习微分方程模型的建立、求解方法、分析结果及解决实际问题的全过程。
1．知道求解微分方程的解析法、数值解法以及图形表示解的方法；
2．理解求微分方程数值解的欧拉方法，了解龙格——库塔方法的思想；
3．熟练掌握使用MATLAB软件的函数求微分方程的解析解、数值解和图形解；
4．通过范例学习怎样建立微分方程模型和分析问题的思想。
§2.1引例
在《大学物理》中，我们曾学习过简谐振动（如：
弹簧振子、单摆）
，那就是一个典型的二阶常微分方程的模型。
这里我们讨论“倒葫芦形状容器壁上的刻度问题”。
对于圆柱形状容器壁上的容积刻度，可以利用圆柱体体积公式：
，其中容器的直径D为常数，体积V与相对于容器底部的任意高度H成正比，因此在容器壁上可以方便地标出容积刻度。
而对于几何形状不规则的容器，比如“倒葫芦形状”的容器壁上如何标出容积刻度呢？
如图所示，建立坐标系，由微元法分析可知：
，其中x表示高度，直径是高度的函数，记为D（x）。
可得微分方程：
如果该方程中的函数D（x）无解析表达式，只给出D（x）的部分测试数据，如何求解此微分方程呢？
h=0.2;
d=[0.04,0.11,0.26,0.56,1.04,1.17];
x
（1）=0;v
（1）=0;
fork=1:
5
x（k+1）=x（k）+h;
v（k+1）=v（k）+（h/2）*（pi/4）*（d（k）^2+d（k+1）^2）;
end
x=x（1:
6）,v=v（1:
6）,
plot（x,v）
x=
Columns1through5
00.20000.40000.60000.8000
Column6
1.0000
v=
Columns1through5
00.00110.00730.03730.1469
Column6
0.3393
§2.2微分方程模型的建立
在工程实际问题中，“改变”、“变化”、“增加”、“减少”等关键词提示我们注意什么量在变化，关键词“速率”、“增长”、“衰变”、“边际的”等常涉及到导数。
我们熟悉的速度公式：
就是一个简单的一阶微分方程。
微分方程是指含有导数或微分的等式。
一般形式：
常用的建立微分方程的方法有：
运用已知物理定律；利用平衡与增长式；运用微元法；应用分析法。
2.2.1运用已知物理定律
建立微分方程模型时应用已知物理定律，可事半功倍。
例2.1一个较热的物体置于室温为180C的房间内，该物体最初的温度是600C，3分钟以后降到500C。
想知道它的温度降到300C需要多少时间？
10分钟以后它的温度是多少？
牛顿冷却（加热）定律：
将温度为T的物体放入处于常温m的介质中时，T的变化速率正比于T与周围介质的温度差。
分析：
假设房间足够大，放入温度较低或较高的物体时，室内温度基本不受影响，即室温分布均衡，保持为m，采用牛顿冷却定律是一个相当好的近似。
建立模型：
设物体在冷却过程中的温度为T（t），t≥0，
根据牛顿加热（冷却）定律：
，建立微分方程
（2.1）
其中参数k>0，m=18。
2.2.2利用平衡与增长式
许多研究对象在数量上常常表现出某种不变的特性，如封闭区域内的能量、货币量等。
利用变量间的平衡与增长特性，可分析和建立有关变量间的相互关系。
此类建模方法的关键是分析并正确描述基本模型的右端，使平衡式成立。
例2.2战斗模型：
两方军队交战，希望为这场战斗建立一个数学模型，应用这个模型达到如下目的：
1.预测哪一方将获胜？
2.估计获胜的一方最后剩下多少士兵？
3.计算失败的一方开始时必须投入多少士兵才能赢得这场战斗？
解：
模型建立：
设x（t）：
t时刻X方存活的士兵数
y（t）：
t时刻Y方存活的士兵数
假设：
1）双方所有士兵不是战死就是活着参加战斗，x（t）与y（t）都是连续变量；
2）Y方军队的一个士兵在单位时间内杀死X方军队a名士兵；
3）X方军队的一个士兵在单位时间内杀死Y方军队b名士兵；
{Δt时间内X军队减少的士兵数}={Δt时间内Y军队消灭对方的士兵数}
即有Δx=－ayΔt
同理Δy=－bxΔt
令
，得到微分方程组：
（2.2）
2.2.3微元法
基本思想：
通过分析研究对象的有关变量在一个很短时间内的变化情况建立微分方程。
例2.3一个高为2m的球体容器里盛了一半的水，水从它的底部小孔流出，小孔的横截面积为1cm2。
试求放空容器所需要的时间。
解：
对孔口的流速做两条假设：
（1）t时刻的流速依赖于此刻容器内水的高度h（t）.
（2）整个放水过程无能量损失。
分析：
放空容器意味着
模型建立：
由流体力学知：
水从孔口流出的流速Q为通过“孔口横截面的水的体积V对时间t的变化率”，即
（孔口流速公式）（2.3）
S—孔口横截面积（单位：
cm2）
h（t）—水面高度（单位：
cm）
t—时间（单位：
s）
当S=1cm2，有
。
2.2.4分析法
基本思想：
根据对现实对象特性的认识，分析其因果关系，找出反映内部机理的规律。
例2.4独家广告模型广告是调整商品销售的强有力的手段，广告与销售量之间有什么内在联系？
如何评价不同时期的广告效果？
解：
1．分析广告的效果，可做如下的条件假设：
①商品的销售速度会因广告而增大，当商品在市场上趋于饱和时，销售速度将趋于一个极限值；
②商品销售率（销售加速度）随商品销售速度的增高而降低。
2．符号说明
A（t）—t时刻的广告费用
S（t）—t时刻商品的销售速度；
M—销售饱和水平，即销售速度的上限；
λ—衰减因子，广告作用随时间的推移而自然衰减的速度，λ＞0；
p—响应系数，表征A（t）对S（t）的影响力。
3．模型建立
选择如下广告策略，t时刻的广告费用为：
建立微分方程：
（2.4）
模型分析：
是否与假设相符？
§2.3微分方程求解方法
2.3.1解析解法
解析解法只能解决一些特殊微分方程，这些方法主要针对：
一阶特殊的微分方程：
如使用分离变量法、方程变换法、线性方程的常数变易法或公式法求解。
二阶或高阶常系数线性微分方程的特征根法。
在高等数学的教程中有专门介绍。
下面着重介绍微分方程的数值解法。
2．3．2数值解法
微分方程的数值解法是解决某些实际问题中经常使用的方法。
设待求解的定解问题为
求该问题数值解法的基本过程如下：
引入自变量取值点序列
，定义
为步长，常用定步长（
与n无关，为常数），其精确解记为
，一般难以得到。
为了寻求
的近似值
，设想根据一定的原理，结合当前得到近似解，近似地表示该点或前一点的导数值，由此推出计算
的迭代公式。
因此数值解法一般只能得到微分方程的近似解
。
下面介绍两个微分方程中最常用的数值解法。
1.欧拉方法
这是一种最简单的解微分方程的数值方法：
就是在小区间[xn,xn+1]上用差商代替微商，可以得到近似的表达式
若f（x,y）中的x取左端点
，结合已经得到的y（xn）的近似值（数值解）yn，即
，有y（xn+1）的近似值为
n=0,1,…
这就是求解微分方程的显式欧拉公式。
也称向前欧拉公式。
向前欧拉法计算简单，易于计算，但精度不高，收敛速度慢
若f（x,y）中的x取右端点
，可得向后欧拉公式如下：
yn+1=yn+hf（xn+1,yn+1），n=0,1,…
称为隐式公式，因为要得出数值解yn+1，就必须求解这个非线性方程，计算比较困难。
如果用将向前和向后欧拉公式加以平均，可得到梯形公式：
该法的计算精度比向前和向后欧拉法都高，但计算和向后欧拉法一样困难。
改进的欧拉算法：
（1）先用向前欧拉法算出yn+1的预测值
，
（2）将预测值
代入梯形公式的右端作为校正，得到yn+1
，n=1,2,…
该式称为改进欧拉公式。
例2.5求解微分方程
y'=-y+x+1,y（0）=1,取步长h=0.1和0.001。
分别用三种数值解法求解，并结合其精确解，对求解误差进行分析比较。
解这是一个一阶线性微分方程，可用解析解法得到其精确解y=x+e-x。
三种数值解如下：
（h=0.1）
1）向前欧拉法：
迭代公式为 yn+1=（1-h）yn+hxn+h，n=0,1,...。
其中y0=y（0）=1。
2）后退欧拉法：
由后退欧拉法隐式公式得yn+1=yn+0.1（-yn+1+xn+1+1），变形为
yn+1=（yn+hxn+1+h）/（1+h）。
3）梯形法：
将隐式梯形公式转化为显示迭代公式如下：
yn+1=（yn+（h/2）*（-yn+xn+xn=1+2））/（1+h/2）。
x1
（1）=0;y1
（1）=1;y2
（1）=1;y3
（1）=1;h=0.1;
fork=1:
10
x1（k+1）=x1（k）+h;
y1（k+1）=（1-h）*y1（k）+h*x1（k）+h;
y2（k+1）=（y2（k）+h*x1（k+1）+h）/（1+h）;
y3（k+1）=（y3（k）+（h/2）*（-y3（k）+x1（k）+x1（k+1）+2））/（1+h/2）;
end
x=0:
0.1:
1;
y=x+exp（-x）;
x1=x1（1:
11）,y=y（1:
11）,y1=y1（1:
11）,y2=y2（1:
11）,y3=y3（1:
11）,
plot（x,y,x1,y1,'k:
',x1,y2,'r--',x1,y3,'g*'）
程序中，x1为自变量，y为精确解，y1、y2、y3分别为向前欧拉法、后退欧拉法和梯形法的解。
结果如下：
x1=
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
y=
Columns1through5
1.00001.00481.01871.04081.0703
Columns6through10
1.10651.14881.19661.24931.3066
Column11
1.3679
y1=
1.0000
1.0000
1.0100
1.0290
1.0561
1.0905
1.1314
1.1783
1.2305
1.2874
1.3487
y2=
1.0000
1.0091
1.0264
1.0513
1.0830
1.1209
1.1645
1.2132
1.2665
1.3241
1.3855
y3=
1.0000
1.0048
1.0186
1.0406
1.0701
1.1063
1.1485
1.1963
1.2490
1.3063
1.3676
图中，蓝色曲线是精确解，红色曲线是向前欧拉法曲线，黑色曲线是向后欧拉法曲线，绿色“*”号为梯形法曲线。
计算结果如下：
表2-1当h=0.1时
精确解
向前欧拉法
后退欧拉法
梯形法
0
1
1
1
1
0.1
1.0048
1
1.0091
1.0048
0.2
1.0187
1.0100
1.0264
1.0186
0.3
1.0408
1.0290
1.0513
1.0406
0.4
1.0703
1.0561
1.0830
1.0701
0.5
1.1065
1.0905
1.1209
1.1063
0.6
1.1488
1.1314
1.1645
1.1485
0.7
1.1966
1.1783
1.2132
1.1963
0.8
1.2493
1.2305
1.2665
1.2490
0.9
1.3066
1.2874
1.3241
1.3063
1
1.3679
1.3487
1.3855
1.3676
当h＝0.001时
x1
（1）=0;y1
（1）=1;y2
（1）=1;y3
（1）=1;h=0.001;
fork=1:
1000
x1（k+1）=x1（k）+h;
y1（k+1）=（1-h）*y1（k）+h*x1（k）+h;
y2（k+1）=（y2（k）+h*x1（k+1）+h）/（1+h）;
y3（k+1）=（y3（k）+（h/2）*（-y3（k）+x1（k）+x1（k+1）+2））/（1+h/2）;
end
x=0:
0.1:
1;
y=x+exp（-x）;
n=1;
fork=1:
11
x1（k）=x1（n）;
y1（k）=y1（n）;
y2（k）=y2（n）;
y3（k）=y3（n）;
n=n+100;
end
x1=x1（1:
11）,y=y（1:
11）,y1=y1（1:
11）,y2=y2（1:
11）,y3=y3（1:
11）,
plot（x,y,x1,y1,'k:
',x1,y2,'r--',x1,y3,'g*'）
x1=
Columns1through5
00.10000.20000.30000.4000
Columns6through10
0.50000.60000.70000.80000.9000
Column11
1.0000
y=
Columns1through5
1.00001.00481.01871.04081.0703
Columns6through10
1.10651.14881.19661.24931.3066
Column11
1.3679
y1=
Columns1through5
1.00001.00481.01861.04071.0702
Columns6through10
1.10641.14861.19641.24911.3064
Column11
1.3677
y2=
Columns1through5
1.00001.00491.01881.04091.0705
Columns6through10
1.10671.14901.19681.24951.3068
Column11
1.3681
y3=
Columns1through5
1.00001.00481.01871.04081.0703
Columns6through10
1.10651.14881.19661.24931.3066
Column11
1.3679
表2-2当h=0.001时
精确解
向前欧拉法
后退欧拉法
梯形法
0
1
1
1
1
0.1
1.0048
1.0048
1.0049
1.0048
0.2
1.0187
1.0186
1.0188
1.0187
0.3
1.0408
1.0407
1.0409
1.0408
0.4
1.0703
1.0702
1.0705
1.0703
0.5
1.1065
1.1064
1.1067
1.1065
0.6
1.1488
1.1486
1.1490
1.1488
0.7
1.1966
1.1964
1.1968
1.1966
0.8
1.2493
1.2491
1.2495
1.2493
0.9
1.3066
1.3064
1.3068
1.3066
1
1.3679
1.3677
1.3681
1.3679
计算结果表明：
当步长h=0.1时，它们的前两位有效数字是精确的；当步长h=0.001时，它们的前四位有效数字是精确的。
说明在迭代中，步长h越小，计算结果越精确。
通过进一步的计算，还可以发现：
迭代离开初始点越远，误差越大。
如上图所示，中间的曲线表示精确解曲线。
2．误差和阶
衡量求解公式好坏的一个主要标准是求解公式的精度，因此引入局部截断误差和阶数的概念。
定义2.1假定
没有误差，即
，用数值方法计算
的误差
，称为该数值方法计算
时的局部截断误差。
为估计欧拉公式的局部截断误差，先将精确解
在
处作泰勒级数展开
（2.12）
对于向前欧拉公式，在
的假定下可记作
（2.13）
由式（2.12）和（2.13）之差可得到
式中h的最低阶项
称为局部截断误差主项，它对于h是2阶的。
同理，对于向后欧拉公式的局部截断误差是
，梯形公式和改进欧拉公式的局部截断误差是
。
求解微分方程的数值计算方法的精度是由局部截断误差中h的阶定义的：
如果一个方法的局部截断误差为
，则称该方法具有p阶精度。
因此向前和向后欧拉方法的精度为1阶，梯形公式和改进欧拉公式的精度为2阶。
3.龙格-库塔方法
数值方法具有越高阶的精度，得到的解就越精确。
我们发现，向前和向后欧拉公式各用了区间
一个端点的导数，具有1阶精度，而梯形和改进欧拉公式用了
两个端点的导数取平均，得到了2阶精度。
因此就引导人们利用
内的若干点的导数，对它们作线性组合得到平均斜率，由此得到更高阶的精度，这就是龙格-库塔方法的基本思路。
Runge-Kutta方法的导出
对于常微分方程的初值问题
的解，在区间
上使用微分中值定理，有
，其中
。
即
设
则
，
即
（2.5）
k可以认为是
在区间上的平均斜率，只要使用不同方法给出
在区间
内平均斜率的近似值k，就可得到不同的数值计算方法。
如：
为向前欧拉方法；
为向后欧拉方法。
1）.二阶R-K方法
一般形式为
（2.15）
若满足
，式（2.15）具有2阶的精度。
如果
时，就得到改进的欧拉方法。
2）.三阶R-K方法
（2.16）
式（2.16）具有3阶的精度。
3）.四阶经典R-K方法
（2.17）
式（2.17）具有4阶精度。
在MATLAB软件中含有数值求解的系统函数，其实现原理就是龙格－库塔方法，MATLAB求解方法将在2.4节中介绍。
§2.4Matlab软件求解
前面已经介绍了微分方程（组）的数值解法。
这些方法计算工作量大，需要借助于计算机实现。
下面简单介绍Matlab在此领域的应用。
2.4.1解析解
用MATLAB命令dsolve（‘eqn1’,’eqn2’,...）求常微分方程（组）的解析解。
其中‘eqni'表示第i个微分方程，Dny表示y的n阶导数,默认的自变量为t。
1.微分方程
例2.6求解一阶微分方程
（1）求通解
输入：
dsolve（'Dy=1+y^2'）
输出：
ans=
tan（t+C1）
（2）求特解
输入：
dsolve（'Dy=1+y^2','y（0）=1','x'）
指定初值为1，自变量为x
输出：
ans=
tan（x+1/4*pi）
例2.7求解二阶微分方程
输入:
dsolve（'D2y+（1/x）*Dy+（1-（1/2）^2/x^2）*y=0','y（pi/2）=2,Dy（pi/2）=-2/pi','x'）
ans=
2^（1/2）*pi^（1/2）/x^（1/2）*sin（x）
化简输出结果，输入：
pretty（ans）
1/21/2
2pisin（x）
-----------------
1/2
x
即：
2．微分方程组
例2.8求解 df/dx=3f+4g; dg/dx=-4f+3g。
（1）通解:
[f,g]=dsolve（'Df=3*f+4*g','Dg=-4*f+3*g'）
f=
exp（3*t）*（C1*sin（4*t）+C2*cos（4*t））
g=
exp（3*t）*（C1*cos（4*t）-C2*sin（4*t））
特解：
[f,g]=dsolve（'Df=3*f+4*g','Dg=-4*f+3*g','f（0）=0,g（0）=1'）
f=
exp（3*t）*sin（4*t）
g=
exp（3*t）*cos（4*t）
2.4.2数值解
在微分方程（组）难以获得解析解的情况下，可以用Matlab方便地求出数值解。
格式为:
[t,y]=ode23（'F',ts,y0,options,p1,p2,...）
注意：
Ø微分方程的形式：
y'=F（t,y），t为自变量，y为因变量（可以是多个，如微分方程组）；
Ø[t,y]为输出矩阵，分别表示自变量和因变量的取值；
ØF代表微分方程组的函数名（m文件，必须返回一个列向量）；
Øts的取法有几种，
（1）ts=[t0,tf]表示自变量的取值范围，
（2）ts=[t0,t1,t2,…,tf],则输出在指定时刻t0,t1,t2,…,tf处给出，（3）ts=t0:
k:
tf,则输出在区间[t0,tf]的等分点给出；
Øy0为初值条件；
Øoptions用于设定误差限（缺省是设定相对误差是10^（-3），绝对误差是10^（-6））；
Øp1,p2,...用于传递附加的参数值。
ode23是微分方程组数值解的低阶方法，ode45为中阶方法，与ode23类似。
比如例2.7的数值解：
解：
令y1=y,y2=y1’,将二阶微分方程转化为一阶微分方程组
y1’=y2
y2’=-y2/x+（（n/x）2-1）y1
首先建立M-文件函数:
functionf=jie3（x,y）
f=[y
（2）;-y
（2）/x+（（1/2）^2/x^2-1）*y
（1）];
计算:
[x,y]=ode23（'jie3',[pi/2,pi],[2,-2/pi]）
x=
1.5708
1.6074
1.7645
1.9215
2.0786
2.2357
2.3928
2.5499
2.7069
2.8640
3.0211
3.1416
y=
2.0000-0.6366
2.早期介入原则；1.9758-0.6869
3）选择价值。
选择价值（OV）又称期权价值。
我们在利用环境资源的时候，并不希望它的功能很快消耗殆尽，也许会设想未来该资源的使用价值会更大。
1.8518-0.8879
1.6982-1.0631
2.环境影响报告表的内容1.5192-1.2108
1.建设项目环境影响评价机构的资质管
升级会员 