打靶法求边值问题.docx

《打靶法求边值问题.docx》由会员分享，可在线阅读，更多相关《打靶法求边值问题.docx（31页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、本科毕业论文设计论文设计题目：打靶法求边值问题 学 院：理学院 专 业：数学应用数学 班 级：091 学 号：0907010228 学生姓名：钟玲声 指导教师：汪萌萌2021年 4月 21日第 28 页打靶法求边值问题目 录摘要：1引言：2第一章 常微分方程初值问题的解法31.1 常微分方程的离散化31.2 欧拉Euler方法41.3 改良的Euler方法61.4 龙格库塔RungeKutta方法71.5 4阶龙格库塔公式91.6 线性多步法9第二章 边值问题的数值解法112.1 打靶法112.2 差分法15第三章 Matlab数值解163.1 常微分方程的解法163.2 打靶法的matlab

2、实现23致 谢：27主要参考文献27摘要常微分方程在很多领域都有非常重要的应用，然而很多常微分方程的解是无法用解析解写出的，因而要借助于数值方法。本文介绍了常微分方程边值问题的常见解法，例如：欧拉法，龙格库塔法等。而对于常微分方程边值问题，常见的解法有打靶法、有限差分法与有限元法等。在本文中，我们重点介绍了打靶法，并给出了相关算法，然后结合实例编写程序进展了上机实验。关键词：常微分方程，初值问题，边值问题，打靶法AbstractOrdinary differential equations play an important role in different areas. However,

3、most equations cannot be expressed analytically. We need to use numerical methods. In this paper, we discuss the methods of solving initial value problem (IVP) , such as Euler method Runge-Kutta method. For boundary value problem (BVP), shooting method, finite difference method (FDM) and finite elem

4、ent method(FEM) are presented. We mainly discuss shooting method and give the algorithm. Numerical experiment is presented in the end of the paper. Keywords: Ordinary differential equations, initial value problem, boundary value problem, the shooting method引言虽然常微分方程理论开展已经有几百年，但目前仍然在开展中。特别是最近三十年，常微分方

5、程迎来了开展的顶峰。常微分方程边值问题是常微分方程理论的重要组成局部, 在众多科学技术领域中有着特别广泛的应用。打靶法是求解常微分方程边值问题的一种数值方法，它的根本思想是将微分方程的边值问题转化为初值问题来求解,它的比拟突出的特点是精度很高,程序很简单,实用性很强。边值问题：对n阶常微分方程 如果能在不同的两点与处，唯一地刻画n个附加条件，并且在区间上求解，那么称此为边值问题。在微分方程中，所谓的边值问题就是我们给定的一个微分方程与一组被我们称之为边界条件的约束条件。边值问题的解一般情况下是符合特定的约束条件的微分方程的解。我们在求解这个微分方程时，除了给出方程的本身，往往还需要提供一定的定

6、解条件。最常见的就是给出初值问题，也就是说给出的定解条件为初始条件；但是也有一些情况，定解条件要求我们考虑所讨论区域的边界，比方说在一个给定区间讨论时，把定解条件在区间的两个端点给出，给定的这种定解条件就被我们称之为边界条件，及之相应的定解问题我们就称之为边值问题。第一章 常微分方程组初值问题的解法1.1 常微分方程的离散化下面主要讨论一阶常微分方程的初值问题，它的一般形式是 (1)在下面的讨论中，总假定函数连续，且满足Lipschitz条件，也就是存在常数，使得那么，根据常微分方程理论知，初值问题(1)的解存在并且唯一.所谓数值解法，就是求问题1的解在假设干点处的近似值的方法，称为问题1的数

7、值解，称为由到的步长.今后如无特别说明，我们总取步长为常量.建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：如果用向前差商代替代入1中的微分方程，那么得化简得如果用的近似值代入上式右端，所得结果作为的近似值，记为，那么有 2这样，问题1的近似解可通过求解下面的问题 3得到，按式3由初值可逐次算出.式子3是个离散化的问题，称为差分方程初值问题.需要说明的是，用不一样的差商近似导数，将得到不一样的计算公式.1.1.2 用数值积分方法将问题1的解表成积分形式，用数值积分方法离散化.例如，对微分方程两端积分，得到 (4)右边的积分用矩形公式或梯形公式计算.1.1.3 Taylor多项式近似

8、将函数在处展开，取一次Taylor多项式近似，那么得再将的近似值代入上式右端，所得结果作为的近似值，得到离散化的计算公式上面的三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式.其中的Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差.1.2欧拉Euler方法1.2.1 Euler方法Euler 方法就是用差分方程初值问题3的解来近似微分方程初值问题1的解，即由公式3依次算出的近似值。这组公式求问题1的数值解称为向前Euler公式.如果在微分方程离散化时，用向后差商代替导数，也就是，那么得计算公式 5用这组公式求问题1的数值解称为向后Euler公式.

9、，因此是隐式公式，一般要用迭代法求解，迭代公式通常为 61.2.2 Euler方法的误差估计 对于向前Euler公式3我们看到，因为时公式右端的都是近似的，所以用它计算的会有累积误差，分析累积误差比拟复杂，这里先讨论比拟简单的所谓局部截断误差.假设用3式时右端的没有误差，即，那么由此算出 7局部截断误差指的是，按7式计算由到这一步的计算值及准确值之差.为了估计它，根据Taylor展开得到的准确值是 87、8两式相减注意到得到 9也就是局部截断误差是阶的，而数值算法的精度定义为： 如果一种算法的局部截断误差为，那么称该算法具有阶精度.显然越大，方法的精度越高。式9说明，向前Euler方法是一阶方

10、法，因此它的精度不高。1.3 改良的Euler方法1.3.1 梯形公式利用数值积分方法将微分方程离散化的时候，如果用梯形公式计算式(4)中之右端积分，即并用代替，那么得计算公式这就是求解初值问题1的梯形公式.直观上容易看出，用梯形公式计算数值积分要比矩形公式好.梯形公式为二阶方法。 梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为 (10) 由于函数关于满足Lipschitz条件，容易看出其中为Lipschitz常数.因此，当时，迭代收敛.但是这样做计算量较大.如果实际计算时精度要求不太高，用公式10求解时，每步可以只迭代一次，由此导出一种新的方法改良Euler法.1.3.2 改良Eule

11、r法按照式5计算问题1的数值解时，如果每步只迭代一次，相当于将Euler公式及梯形公式结合使用：先用Euler公式求的一个初步近似值，称为预测值，然后用梯形公式校正求得近似值，即 11式11称为由Euler公式与梯形公式得到的预测校正系统，也叫改良Euler法.为了便于编制程序上机，式子11常改写成 12改良Euler法是二阶方法.1.4 龙格库塔RungeKutta方法回到Euler方法的根本思想用差商代替导数.实际上，按照微分中值定理应有注意到方程就有 13不妨记，称为区间，13式就对应地导出一种算法.向前Euler公式简单地取为取，的平均值，其中，这种处理提高了精度.如上分析启示我们，在

12、区间内多取几个点，将它们的斜率加权平均作为库塔方法的根本思想. 首先不妨在区间内仍取2个点，仿照13式用以下形式试一下 14其中，因为，所以14可以化为 15其中在点作了Taylor展开.15式又可表为注意到中，可见为使误差，只须令 ， 16待定系数满足16的15式称为2阶龙格库塔公式.由于16式有4个未知数而只有3个方程，所以解不唯一.不难发现，假设令 ，即为改良的Euler公式.可以证明，在内只取2点的龙格库塔公式精度最高为2阶.1.5 四阶龙格库塔公式要进一步提高精度，必须取更多的点，如取4点构造如下形式的公式： 17其中待定系数共13个，经过及推导2阶龙格库塔公式类似、但更复杂的计算，

13、得到使局部误差的11个方程.取既满足这些方程、又较简单的一组，可得 18这就是常用的4阶龙格库塔方法简称RK方法.1.6 线性多步法以上所介绍的各种数值解法都是单步法，这是因为它们在计算时，都只用到前一步的值，单步法的一般形式是 19其中称为增量函数，例如Euler方法的增量函数为，改良Euler法的增量函数为怎么样通过较多地利用前面的信息，如，来构造高精度的算法计算，这就是多步法的根本思想.经常使用的是线性多步法.让我们试验一下，即利用计算的情况.从用数值积分方法离散化方程的4式出发，记，式中被积函数用二节点，的插值公式得到因，所以是外插. 20此式在区间上积分可得于是得到 21注意到插值公

14、式20的误差项含因子，在区间上积分后得出，故公式21的局部截断误差为，精度比向前Euler公式提高1阶.假设取可以用类似的方法推导公式，如对于有 22其局部截断误差为.如果将上面代替被积函数，取时得到的是梯形公式，取时可以得到 23及22式相比，虽然其局部截断误差仍为，但因它的各项系数绝对值大为减小，误差还是小了.当然，23式右端的未知，需要如同向后Euler公式一样，用迭代或校正的方法处理.第二章 边值问题的数值解法2.1 打靶法解两点边值问题 1的打靶法实质上是把边值问题化为初值问题来解，我们设法确定的值，使得初值问题 2的解在的值满足或者是其中为我们允许的误差界限. 这样，我们把看作边值

15、问题1的近似解.为此，可以采用逐次逼近的方法来实现.假设为边值问题的解，我们估计的值为之后，求解初值问题在2中令.这样得到的解为，并且计算得到.一般地，.但是如果或，那么把作为边值问题1的近似解；否那么的话，那么做适当的调整，例如取，那么在2中令，假设或，那么作为边值问题1的近似解；否那么，再做适当的调整.如此重复计算，直至或时，便以作为边值问题1的近似解，参考下列图：上面的积分曲线的过大，下面的曲线的.当时，从初值问题计算得到解.自然地，我们希望因此，确定的问题可以归纳为求方程 3的近似根.3是一个非线性方程.二分法，割线法与牛顿法都可以用来求解方程3.假设用弦割法来求解方程3，我们需要选取

16、初始近似，由公式 4生成序列.按照4的方式求，直到为止，其中为允许的误差界.令.给定初值，误差容限，最大迭代次数.取，解初值问题2得解在的近似值：如果，那么输出作为初值问题1的解在的近似值；停机令，解初值问题2得解在的近似值：如果，那么输出作为初值问题1的解在的近似值；停机.对重复做。由4公式计算；令，解初值问题2的解在的近似值：如果，那么输出作为初值问题1的解在的近似解；停机. 5 6 7打靶法的根本过程如下：选取参数，令.设，那么可以由6式确定，从而得到初值问题 8求初值问题8的解代入7式左端得如果，那么为所求解的边值问题的近似解；如果，那么令设以取代解初值问题8，求得解后在代入7式左端得

17、从出发，由割线法 9产生，以取代解决初值问题8确定，如此继续进展下去，求得序列，直到为止，为误差容限.如果，那么5，6，7为第二边值问题，处理方法完全类似.2.2 差分法在这里我们介绍一种求解1的简单方法，我们以最简单的二阶常微分方程为例来说明边值问题的数值解法.解法可以用数值微分公式替代导数，将其变成代数方程，然后我们再求解，这种求解的方法我们通常把它叫做差分法. 即中心差商 即二阶差商令并略去与用差商表示中的二次项，然后代入1，这样我们就得到： 10其中，.从10中解出，即为所求的常微分方程的计算解.一般情况下，10是一个非线性方程组，求解也比拟困难.但是，如果微分方程1中第一式的右端为的

18、线性函数时，那么方程组10就可以转化成线性方程组，求解问题也就简单多了.第三章 Matlab数值解3.1 常微分方程的解法在Matlab的工具箱中，给我们提供了很多求解常微分方程的功能函数.如ode45，ode23，ode113，其中ode45采用四五阶RK方法，是解非刚性常微分方程的首选方法，ode23采用二三阶RK方法，ode113采用的是多步法，效率一般比ode45高.Matlab的工具箱中没有Euler方法的功能函数.对简单的一阶方程的初值问题改良的Euler公式为一阶微分方程初值问题的求解在Matlab中非常简单明了并且容易求解的方法，而且由于Matlab中内嵌的微分方程初值问题的求

19、解方法是变步长的龙格库塔法，故而求解的结果是值得我们相信的.对于微分方程，直接利用Matlab语言的内嵌函数ode45，按照调用格式x,y=ode45(,指函数名或所建立的M文件的文件名，为离散节点向量,为初值.例：求解一阶微分方程初值问题： 1先编写M文件fun.m:Function f=fun(x,y)f=y-2*x/y;命令窗口写入：x=0:0.2:1;x,y=ode45(fun,x,1)求解结果为：这个结果及四阶龙格库塔法求解结果一样，比改良欧拉法求解结果要准确得多.运用内嵌函数直接求解最大的优点是简洁方便，只要依照函数的调用格式调用即可，但是对于理解函数值解法的思想没有太大的帮助，故

20、而要表达数值解法的具体思想，可以通过编写程序来实现.3.1.2 根据数值解法的思想编程实现对于一阶微分方程的初值问题 2运用改良欧拉法与四阶龙格库塔法求解程序如下：(1) 改良欧拉法建立如下M文件：a=0;b=5;n=120;xn=0;yn=2;h=(b-a)/n;for i=1:nyp=yn-h*xn*yn*yn;xn=xn+h;yc=yn-h*(xn*yp*yp);yn=(yp+yc)/2;end迭代次数n=120时求解结果：yn=0.450;假设n=5000时,求解结果yn=0.101.标准四阶龙格库塔法对于上述初值问题建立如下M文件：a=0;b=5;n=20;xn=0;yn=2;h=(

21、b-a)/h;for i=1:n； k1=-h*xn*yn*yn; xn=xn+h/2; k2=-h*xn*(yn+k1/2)*(yn+k1/2); k3=-h*xn*(yn+k2/2)*(yn+k2/2); xn=xn+h/2; k4=-h*xn*(yn+k3/2)*(yn+k3); yn=yn+(k1+2*k2+2*k3+k4)/6;end当迭代次数n=120的时候，获得计算结果：yn=0.821，比运用改良欧拉法迭代120步的效果都要好的多；儿当迭代次数n=120时，迭代结果yn=0.295根本已经到达准确解.3.1.3 结果分析对问题1的求解简单方便，但难于表达求解的具体思想：对于问题

22、2的求解运用变成实现，从实验角度帮助理解理论知识，不仅可以深入掌握求解方法的思想，而且还可以提高实践动手能力.，改良欧拉法在迭代到120步的时候误差还较大，迭代到5000步的时候误差较小，但是及准确解的误差仍然有出入，而且从迭代过程来看，迭代次数n越大，收敛速度越慢，效果不明显.而四阶龙格库塔法迭代20步的结果远远好于改良欧拉法120步的结果，当迭代到120步时已经跟准确解非常接近.故而龙格库塔法迭代不仅得到的解比改良欧拉法准确得多，而且收敛速度非常快，确实是一种很好的算法.但是唯一缺乏的是计算量大，这一矛盾很难解决.3.2 打靶法的Matlab实现 我们考虑两个边界条件的二阶常微分方程其中是

23、给定的常数，是关于的未知的函数，是一个函数的微分方程的指定者，实现了及值下.例如，与.因为两个独立的边界条件，上述两点边值问题是比拟难解决的问题.“打靶法是一个初始值问题代替上述问题的根本思想.当然，我们并不知道在时的导数值，但是我们可以猜想，然后进一步完善猜想的迭代.更确切的说是，我们对为未知，并使用割线法或者牛顿法或者其他的方法求解非线性方程组来确定.我们使用与表示相对于是准确的，如果是正确的选择.通过以下定义：当给定的值，然后就可以求解上述关于及是不等的，由于.但是如果等于，那么我们认为=.因为我们不知道，我们确定它在边界条件.即，我们要求解使得：如果我们找到了使得，这意味着.因此，在与

24、，满足两个一样的边界条件，换句话说，.因此，的解必须满足.如果我们可以写出初始值问题任意的解析式，我们可以写一个公式.当然，这一般是不可能的.然而，没有一个解析公式，我们仍然可以求解，初始值问题的数值方法可以用来找到一个近似解因此.最简单的方法是使用割线法.为了到达这个目的，我们需要两个初始猜想：与.我们也可以使用牛顿法：我们需要一种方法来计算导数，由于，我们有如果我们定义了，我们有一下的初始值问题：在这里，对于，与是系统的微分方程.现在，我们设.这是打靶法的方法，其中包括牛顿法求解的算法：=初始猜想值.for j=0,1,2,. 求解如下的微分方程系统，求解区间从x=a到x=b =f(x,u

25、,) = =0 =1.令如果我们想使用在前面解决了的上述系统中的两个第二阶方程与的方法，我们需要引入一个向量与微分方程为的向量.初始条件是.打靶法也使用于特征值问题：在满足条件或者跟普遍，在是均匀的，即，对于任何常数.注意总是一个上述边值问题的解.事实上，一个特征值问题是一个特殊的边界值问题满足始终是一个解决方案，有一个称为的方程参数或边界条件.特征值问题是确定非零解为，其中是特征值问题的特征.通常，有许多可能是无限的的特征值与本征函数.利用打靶法，我们考虑初始值问题其中，我们用符号，但与代表第一与第二阶导数取决于.对于任何给定的满足条件那么是我们正在寻找的特征函数，是相应的特征值.因此，我们

26、只需要使用割线法或者牛顿法求解的方程.如果使用割线法，我们只需要解决不同的迭代器的初值问题.如果用牛顿法，我们必须得到关于的.因此，我们需要我们需要一个初始值问题的.这可以通过偏导数相对于的出数值问题，得到.我们有请注意，我们已经使用链式法那么微积分为得到方程.现在，你可以解决初始值问题与一起，然后评估为任何给定的.3.2.2 实例(1) 用打靶法算法编写程序上机求解下面的实例：解 编写算法如下(我们用D2u表示u的2阶导数，用Du表示u的一阶导数；用D2v表示v的2阶导数，用Dv表示v的一阶导数)：t0=0.for j=0,1,2,. solve the following system n

27、umerically from x=-1 to x=1 D2u=f(x,u,Du) u(x=-1)=1 Du(x=-1)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=-1)=0 Dv(x=-1)=1.set t(j+1)=t(j)-(u(x=1)-1)/v(x=1).plot(t,y)用上面的程序上机进展测试，就会得到一个很接近准确解的一个解.(2) 用打靶法算法编写程序上机求解下面的实例：解 编写算法如下(我们用D2u表示u的2阶导数，用Du表示u的一阶导数；用D2v表示v的2阶导数，用Dv表示v的一阶导数)：t0=0.for j=0,1,2,. solv

28、e the following system numerically from x=0 to x=/4 D2u=f(x,u,Du) u(x=0)=0 Du(x=0)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=0)=0 Dv(x=0)=1.set t(j+1)=t(j)-(u(x=/4)-0)/v(x=/4).plot(t,y)(3) 用打靶法算法编写程序上机求解下面的实例：解 编写算法如下(我们用D2u表示u的2阶导数，用Du表示u的一阶导数；用D2v表示v的2阶导数，用Dv表示v的一阶导数)：t0=0.for j=0,1,2,. solve the fo

29、llowing system numerically from x=0 to x=1 D2u=f(x,u,Du) u(x=0)=0 Du(x=0)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=0)=0 Dv(x=0)=1.set t(j+1)=t(j)-(u(x=1)-0)/v(x=1).plot(t,y)(4) 用打靶法算法编写程序上机求解下面的实例：解 编写算法如下(我们用D2u表示u的2阶导数，用Du表示u的一阶导数；用D2v表示v的2阶导数，用Dv表示v的一阶导数)：t0=0.for j=0,1,2,. solve the following sys

30、tem numerically from x=0 to x=1 D2u=f(x,u,Du) u(x=0)=0 Du(x=0)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=0)=0 Dv(x=0)=1.set t(j+1)=t(j)-(u(x=1)-1)/v(x=1).plot(t,y)(5) 用打靶法算法编写程序上机求解下面的实例：解 编写算法如下(我们用D2u表示u的2阶导数，用Du表示u的一阶导数；用D2v表示v的2阶导数，用Dv表示v的一阶导数)：t0=0.for j=0,1,2,. solve the following system numeric

31、ally from x=0 to x=1 D2u=f(x,u,Du) u(x=0)=1 Du(x=0)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=1)=0 Dv(x=1)=1.set t(j+1)=t(j)-(u(x=1)-2)/v(x=1).plot(t,y)(6) 用打靶法算法编写程序上机求解下面的实例：解 编写算法如下(我们用D2u表示u的2阶导数，用Du表示u的一阶导数；用D2v表示v的2阶导数，用Dv表示v的一阶导数)：t0=0.for j=0,1,2,. solve the following system numerically from x

32、=1 to x=2 D2u=f(x,u,Du) u(x=1)=0 Du(x=1)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=1)=0 Dv(x=1)=1.set t(j+1)=t(j)-(u(x=2)-2)/v(x=2).plot(t,y)(7) 用打靶法算法编写程序上机求解下面的实例：t0=0.for j=0,1,2,. solve the following system numerically from x=1 to x=2 D2u=f(x,u,Du) Du(x=1)=t(j) D2v= fu(x,u,Du)v+fDu(x,u,Du)Dv v(x=1

33、)=0 Dv(x=1)=1.set t(j+1)=t(j)-(u(x=2)-ln2)/v(x=2).plot(t,y)致 谢：在本论文完成之际，我要由衷感谢汪萌萌教师在课题设计与论文写作上的精心指导，同时对所有帮助过我的教师以及同学致以最真诚的感谢。主要参考文献：1 施吉林，刘淑珍，陈桂芝.计算机数值方法(M).北京.高等教育出版社.2021.165.2 安乐.打靶法在常微分方程边值问题中的一些应用J.科技广场.2021年.05期.3 夏必腊;王金山.自共轭常微分方程边值问题及变分问题的等价性J.大学数学.2021(03)4 王国栋.常微分方程中数学建模思想的教学探讨J.中国科教创新导刊.2021年.22期.5 王国栋.常微分方程中数学思想方法教学探讨J.科教新报(教育科研).2021年.22期.6 朱美玲.常微分方程在数学建模中的应用J.太原城市职业技术学院学报.2021.05期.7 吴德林.常微分方程教学探索J.读及写(教育教学刊).2021年04期.8 张建龙.常微分方程教学改革的探索J.才智.2021.16期.9 李培峦.几类非线性微分方程边值问题解的存在性及多解性研究D.中南大学.2021年.10 连海强.几类非线性微分方程边值问题的解及应用D.曲阜师范大学.2021年.