第十一章 常微分方程边值问题的数值解法

1、第十一章 常微分方程边值问题的数值解法工程技术与科学实验中提出的大量问题是常微分方程边值问题本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法.11.1 引言在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程, , (11.1.1) 在如下三种边界条件下的定解问题:第一种边界条件: , (11.1.2) 第二种边界条件: , (11.1.2)第三种边界条件: , (11.1.13)其中. 常微分方程边值问题有很多不同解法, 本书仅介绍打靶方法和有限差分方法.11.2 打靶法对于二阶非线

2、性边值问题 (11.2.1)打靶法近似于使用初值求解的情况 我们需要利用一个如下形式问题初值解的序列：， （11.2.2）引进参数以近似原边界值问题的解选择参数，以使：， （11.2.3）其中定义为初值问题（11.2.2）在时的解，同时定义为边值问题（11.2.1）的解首先定义参数，沿着如下初值问题解的曲线，可以求出点对应的初始正视图 （11.2.4）如果不严格收敛于，那么我们选择等值以修正近似值，直到严格逼近为了取得合适的参数，现在假定边值问题（11.2.1）有唯一解，如果定义为初始问题（11.2.2）的解，那么可由下式确定： （11.2.5）由于这是一个非线性方程，我们可以利用Newton

3、法求解首先选择初始值，然后由下式生成序列，此处， （11.2.6）同时要求求得，因为的表达式未知，所以求解这个有一点难度；我们只能得到这么一系列的值假如我们如下改写初值问题（11.2.2），使其强调解对和的依赖性，（11.2.7）保留初始记号以显式与的微分相关既然要求当时的值，那么我们需要求出表达式（11.2.7）关于的偏导数过程如下： 又因为跟相互独立，所以当上式如下；（11.2.8）初始条件为：如果简单地用定义，并且假定和的微分阶翻转，（11.2.8）转化为初值问题：，（11.2.9）因此，对于（11.2.2）和（11.2.9）式的每次迭代需要求解两个初值问题那么从（11.2.6）式可得：

4、， （11.2.10）事实上，这些初值问题很难精确求解，而将这些解近似为一个初值问题的解同样，我们可以按以上步骤考虑对于三阶非线性边值问题的打靶法算法对于三阶非线性边值问题(11.2.11)转变形式：,(11.2.12)选择参数，以使：， （11.2.13）其中定义为初值问题（11.2.12）在时的解，同时定义为边值问题（11.2.11）的解定义参数，沿着如下初值问题解的曲线，可以求出点对应的初始正视图（11.2.14）如果不严格收敛于，那么我们选择等值以修正近似值，直到严格逼近为了取得合适的参数，现在假定边值问题（11.2.11）有唯一解，如果定义为初始问题（11.2.12）的解，那么可由下

5、式确定： （11.2.15）利用Newton法求解首先选择初始值，然后由下式生成序列，此处， （11.2.16）同时要求求得，因为的表达式未知，所以只能得到这么一系列的值改写初值问题（11.2.12），使其强调解对和的依赖性 （11.2.17）要求当时的值，那么我们需要求出表达式（11.2.17）关于的偏导数所以当上式如下； （11.2.18）初始条件为：对于第三种边界条件, 注意到 , 取 得到, , (11.2.19) 设它的解为, 代入第二个边界条件得到即 当 时, 即为所求的解. 这样,同第一种边界条件一样, 可以利用打靶方法求解.例11.1 利用打靶法求解两点边值问题解：(1) 为应

6、用打靶方法，需要假定初值来先求解初值问题，取初始参数 ， 解 令，则上述初值问题变为， 由标准的R-K 方法，取，可得。（2）再令， 解， 取。仍由标准的R-K 方法，可得。（3）令再解， 取，得。（4）令，得（5）类似有，因为， 可以认为即为所求。11.3 有限差分方法有限差分方法常用来求解常微分方程边值问题, 下面就线性微分方程为例进行讨论. 设微分方程具有形式:, , (11.3.1) , . (11.3.2)不失一般性, 将求解区间等份, 分点为, , , 步长, 这里 称为结点.在内部结点上, 记 , , , , 并分别用一阶和二阶中心差分来代替一阶导数和二阶导数, ,.将上面各式代

7、入方程(11.3.1), 略去项, 得到关系式:,将上式改写成, , (11.3.3)在上式中, , , (11.3.4) 利用边界条件 , 可以得到如下形式的方程组: (11.3.5)其中 (11.3.6) (11.3.7) 系数矩阵(11.3.6)为n-1阶三对角阵, 当系数矩阵非奇异时,方程组(11.3.5) 有唯一解,可以利用追赶法求解.对于第二种边界条件: , (11.3.8)和第三种边界条件: , (11.3.9)其中. 由于边界条件中包含了导数, 故边界条件也必须用差商近似来表示. 因为不能利用区间以外的点,所以在引进(即 或 (即 ) 的近似式时,就不能利用中心差商. 如果只要

8、求误差是, 则可以用最简单的近似公式, (11.3.10)要使误差达到, 可以利用Newton等距插值公式,得到如下近似公式:, (11.3.11), (11.3.12)将这些近似式代入边值条件中, 再和方程(11.3.4)联立, 就可以得到对应的差分方程组. 例11.2 利用差分法求解两点边值问题解：将方程离散后写成矩阵形式取， 可以得到三对角方程，下表11-1给出了计算结果。表 11-1取时计算结果，近似值近似值精确值0.11074610744107430.21170111696116950.31287712871128690.41429414286142840.5159761596815

9、9650.61795817950179470.72028520276202740.82301323006230040.926219262152621411.4 应用实例与Matlab11.4.1 MATLAB 关于常微分方程边值问题数值解法的命令常微分方程初、边值问题的符号解法 命令形式1：dsolve(，eqution，) 功能：求常微分方程eqution的解。 命令形式2：dsolve(，eqution，condl，cond2，var，) 功能：求常微分方程eqution的满足初始条件的特解。其中eqution是求解的微分方程或微分方程组，condl，cond2是初始条件或边值，var是自

10、变量.例11.3求解两点边值问题：.解 Matlab命令为 y=dsolve(xD2y-3Dy=x2，y(1)=0，y(5)=0,x)执行结果： y 13x3+125468+31468x4 11.4.2 应用实例当流体以均匀的流速纵掠一平壁时，由于流体粘性的影响，在靠近壁面邻近会形成很大的速度梯度，这就是速度边界层用来描述速度边界层的动量方程经过相似变换可以化为如下常微分方程两点边界值问题(著名的布拉修斯边界层方程) (11.3.13)边界条件为 (11.3.14)这里表示壁面喷注(或吸入). 根据给定的喷注速度值，(常数)表示吸入， (常数)表示喷注.为了利用打靶法求解该问题的解，假设（这里

11、是一个待定常数，需要在计算中估计). 表11-1 给出了数值计算方法。当时取、和.这里是一个常数(估计值)，这样就能根据式 (11.3.13) - (11.3.14), 求得.接着，令具有某一增量，分别计算、和，继续进行下去一直到边界外边缘(在平壁的情况下它相当于). 此时，应该得到. 如果这个要求不能满足，就要修改原来的估计值，直至使这个边界条件得到满足为止.，表11-2 给出了数值计算方法。表11-2 布拉修斯方程式的数值计算方法0N0式(1.31)0.1N+0.1*00+0.1*+0.1*式(1.31)0.2N+0.1*0+0.1*(0+0.1*)0+0.1*+0.1*+0.1*)+0.

12、1*+0.1*式(1.31)50.991.0？在程序中为了能够方便的变换系数，我们把方程中系数设为epsilo，从而把方程改写为： (11.3.13) Matlab 命令： ?itBlasius(0,0,0,0.5,0.00001,0.1,5,0.001)beta=0.435860得出每一步的迭代结果：Matlab 命令：?blasius(0,0,0.435860,0.1,5,0.001,1,1,1,1,1)beta=0.435860表11-2给出了当情况在区间0，5的部分计算结果.Matlab程序function itaconv=Blasius(n,m,beta,dita,maxita,to

13、l,options)format long;%(f+epsilo*f*f=0)%书上是epsilo=0.5%epsilo=1;epsilo=0.5;%if options(5)=1 fprintf(nbeta=%fnttitattfttfttfttfnn,beta);endfound=0;x0=n;y0=m;z0=beta;w0=-x0*z0*epsilo;for ita=0 : dita : maxita switch options(1) case 1 multiplier=0.1;%推荐 otherwise fprintf(n出错了！nn); break end switch optio

14、ns(2) case 1 x=x0+multiplier*y0;%推荐 otherwise fprintf(n出错了！nn); break end switch options(3) case 1 y=y0+multiplier*z0;%推荐 otherwise fprintf(n出错了！nn); break end switch options(4) case 1 z=z0+multiplier*w0;%推荐 otherwise fprintf(n出错了！nn); break end w=-x*z*epsilo; % if options(5)=1 fprintf(t%2.5ft%2.5ft

15、%2.5ft%2.5ft%2.5fn,ita,x0,y0,z0,w0); end% if abs(1-y)-1 fprintf(n满足白拉修斯收敛条件的beta已经找到，beta=%fnn,beta); found=found+1; break%若需找到所有beta, 请注释掉上一行的break endendif found=0 fprintf(n在给定的范围 (%f.%f) 内无法找到满足白拉修斯收敛条件的beta.nn,minbeta,maxbeta);end下面表11-3 给出了当情况在区间0，5的部分计算结果表11-3：当部分计算结果. 0000.332060.20.006640.06

16、6410.331990.40.026560.132770.331470.60.059740.198940.330080.80.106110.264710.3273910.165570.329790.323011.20.237950.395780.316591.40.322980.456270.307871.60.420320.516760.296671.80.529520.574770.2829320.650030.629770.266752.20.78120.681320.248352.40.92230.728990.228092.61.072520.772460.206462.81.230

17、990.811520.1840131.396820.846050.161363.21.569110.876090.139133.41.746960.901770.117883.61.929540.923330.098093.82.116050.941120.0801342.305760.955520.064244.22.498060.966960.050524.42.692380.975870.038974.62.888260.982690.029484.83.085340.987790.0218753.283290.991550.01591小结本章研究了常微分方程边值问题的数值求解方法. 主

18、要介绍了三种边界条件下的定解问题, 通常主要有两大类求解方法, 即:打靶法算法和有限差分方法. 打靶算法对于线性或非线性边界值问题均适应. 边值问题的解可以视为一个 “弹道曲线”, 打靶算法是通过引进初值参数以近似原边界值问题, 即将原微分方程先看成一个初值问题求解, 需要利用一系列初值形式问题的解的序列来逼近, 沿着初值问题解的曲线积分，使该初值问题的解满足另一端点的边界条件, 每一步计算都要进行判断, 通过迭代得到新的初值参数, 以达到满足右端边界条件. 对于不同的边界条件, 需要引入不同形式的初值参数及建立不同代格式. 有限差分方法通常仅适应于求解常微分方程边值问题, 对于求解常微分方程初值问题有时则较困难. 通常通过将求解区间等份, 并分别用一阶和二阶中心差分来代替一阶导数和二阶