常微分方程边值问题的数值解法

1、常微分方程论文题目： 常微分方程边值问题的数值解法 组 长： 数学132文洲 组 员： 数学131王琦 数学132姚瑶 信息132郭斌 院 （系）： 理学院 指导教师： 岳宗敏 时 间： 2015年6月9日 常微分方程边值问题的数值解法摘要：作为一类定解问题，补充条件由以自变量取某些值时，未知函数及其导数的值而定，称其为边值条件。许多物理和数学问题都归结为边值问题。本文介绍边值问题的待定常数法和格林函数。关键词：边值问题 待定常数法 格林函数Abstract: as a kind of definite solution problems, the supplementary conditio

2、ns by took the certain values in the independent variable, the value of the unknown function and its derivative, referred to as boundary value conditions. Many physical and mathematical problems boil down to boundary value problems. In this paper, the boundary value problem of the method of undeterm

3、ined constants and green's function.Keywords: boundary value problem Method of undetermined constants Green's function11.1 引言在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程, , (11.1.1) 在如下三种边界条件下的定解问题:第一种边界条件: , (11.1.2) 第二种边界条件: , (11.1.2)第三种边界条件: , (11.1.13)其中. 常微分方程边值问题有很多不同解法, 这里只介绍打靶方法和有限差分方法.1

4、1.2 打靶算法 将边值问题转化成初值问题来求解，即根据边界条件(11.1.2)，也就是说，反复是调整初始时刻的斜率值，使得初值问题的积分曲线能“命中”。例11.1 利用打靶法求解两点边值问题解：(1) 为应用打靶方法，需要假定初值来先求解初值问题，取初始参数 ， 解 令，则上述初值问题变为， 由标准的R-K 方法，取，可得。（2）再令， 解， 取。仍由标准的R-K 方法，可得。（3）令再解， 取，得。（4）令，得（5）类似有，因为， 可以认为即为所求。11.3 有限差分方法原理：微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替， 这些离散点称作网格的节点；把连续

5、定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组  ， 解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。有限差分方法常用来求解常微分方程边值问题, 下面就线性微分方程为例进行讨论. 设微分方程具有形式:, , (11.3.1) , . (11.3.2)不失一般性, 将求解区间等份, 分点为, , , 步长, 这里 称为结点.在内部结点上, 记 , , , , 并分别用一阶和二阶中心

6、差分来代替一阶导数和二阶导数, ,.将上面各式代入方程(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)其中. 由于边界条件中包含了导数, 故边界条件也必须用差商近似来表示. 因为不能利用区间以外的点,所以在引进(即 或 (即

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

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

9、ondl，cond2是初始条件或边值，var是自变量.例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)这里表示壁面喷注(或吸入). 根据给定的喷注速度值，(常数)表示吸入， (常数)

10、表示喷注.为了利用打靶法求解该问题的解，假设（这里是一个待定常数，需要在计算中估计). 表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+

11、0.1*)0+0.1*+0.1*+0.1*)+0.1*+0.1*式(1.31)50.991.0在程序中为了能够方便的变换系数，我们把方程中系数设为epsilo，从而把方程改写为： (11.3.13) function itBlasius(n,m,minbeta,maxbeta,dbeta,dita,maxita,tol)found=0;fprintf('n白拉修斯方程数值计算迭代程序nn正在计算.nn');for beta=minbeta:dbeta:maxbeta ita=Blasius(n,m,beta,dita,maxita,tol,1,1,1,1,0); if ita&

12、gt;-1 fprintf('n满足白拉修斯收敛条件的beta已经找到，beta=%fnn',beta); found=found+1; break endendif found=0 fprintf('n在给定的范围 (%f.%f) 内无法找到满足白拉修斯收敛条件的beta.nn',minbeta,maxbeta);endMatlab 命令：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

13、,1,1)beta=0.435860表11-2给出了当情况在区间0，4的部分计算结果.Matlab程序function itaconv=Blasius(n,m,beta,dita,maxita,tol,options)format long;epsilo=0.5;if options(5)=1 fprintf('nbeta=%fnttitattfttf''ttf''''ttf''''''nn',beta);endfound=0;x0=n;y0=m;z0=beta;w0=-x0*z0*

14、epsilo;for ita=0 : dita : maxita switch options(1) case 1 multiplier=0.1; otherwise fprintf('n出错了！nn'); break end switch options(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');

15、 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%2.5ft%2.5ft%2.5fn',ita,x0,y0,z0,w0); end if abs(1-y)<tol found=1;itaconv=ita; break end x0=x;y0=y;z0=z;w0=w;endif found itaconv=

16、-1;end下面表11-3 给出了当情况在区间0，5的部分计算结果表11-3：当部分计算结果. 0000.332060.20.006640.066410.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