第十一章 常微分方程边值问题的数值解法汇总
常微分方程的边值问题
![常微分方程的边值问题](https://img.taocdn.com/s3/m/52635828a88271fe910ef12d2af90242a895abfb.png)
常微分方程的边值问题常微分方程的边值问题是指在一定的边界条件下,求解常微分方程的解。
对于常微分方程而言,通常有两种类型的问题:初值问题和边值问题。
初值问题是在给定初始条件下,求解方程的解;而边值问题则是在给定边界条件下,求解方程的解。
边值问题在实际问题中具有广泛的应用,例如求解杆的挠度、求解电路中的电流分布等等。
在这篇文章中,我们将重点讨论边值问题的求解方法及其应用。
首先,我们来回顾一下常微分方程的一般形式:$$%frac{d^2y}{dx^2} = f(x, y, %frac{dy}{dx})$$其中,$f(x, y, %frac{dy}{dx})$是已知函数。
对于一般的边值问题,我们需要给定方程的边界条件,即在一定的边界点上,已知方程的解或者导数值。
通常,边界条件可以分为两类:Dirichlet条件和Neumann条件。
Dirichlet条件是指在边界点上给定方程的解,即$y(a)= %alpha$和$y(b) = ¾ta$,其中$a$和$b$是区间的端点,$%alpha$和$¾ta$是已知的常数。
Neumann条件是指在边界点上给定方程的导数值,即$%frac{dy}{dx}(a) = %alpha$和$%frac{dy}{dx}(b) = ¾ta$。
接下来,我们将介绍边值问题的求解方法。
对于边值问题的求解,最常用的方法是有限差分法。
有限差分法将区间离散化,并用差分代替微分,将微分方程转化为差分方程。
然后,通过求解差分方程,得到方程的数值解。
在有限差分法中,我们首先将区间$[a, b]$离散化为$n$个小区间,即$a = x_0 < x_1 < x_2 < … < x_n = b$,其中$n$是离散点的个数。
然后,我们用$y_i$来表示$y(x_i)$,$y_i’$来表示$%frac{dy}{dx}(x_i)$,并使用差分公式来逼近微分方程中的导数。
常微分方程的边值问题
![常微分方程的边值问题](https://img.taocdn.com/s3/m/906e73501fb91a37f111f18583d049649b660e39.png)
常微分方程的边值问题一、引言在数学中,微分方程是研究自然界中变化和发展的重要工具。
它描述了物体在不同变化条件下的行为规律,并被广泛应用于物理、工程、经济等领域。
边值问题是微分方程中的一个重要分支,它关注的是在一定边界条件下的解。
二、常微分方程常微分方程是指只含有关于一个自变量的一阶或高阶导数的方程。
一般形式为:[F(x, y, y’, y’’, , y^{(n)}) = 0]其中,x是自变量,y是未知函数。
常微分方程的求解可以分为两种类型:初值问题和边值问题。
三、边值问题的定义边值问题是指在一定边界条件下,求解微分方程的解。
对于二阶常微分方程,边值问题的一般形式为:[y’‘(x) = f(x, y, y’), a < x < b, y(a) = , y(b) = ]其中,a和b是给定的边界点,()和()是给定的边界值。
四、边值问题的求解方法边值问题的求解可以分为两种方法:迭代方法和直接方法。
4.1 迭代方法迭代方法是通过不断迭代逼近的方式求解边值问题。
常用的迭代方法有有限差分法和有限元法。
4.1.1 有限差分法有限差分法是一种将微分方程转化为差分方程进行求解的方法。
它将求解域离散化,并通过差分近似来近似微分项,最终通过迭代逼近求得边界值。
有限差分法的基本思想是将求解域划分为若干个离散的网格点,然后使用近似公式将微分项替换为差分项,从而得到差分方程。
通过迭代求解差分方程,最终得到边界条件下的解。
4.1.2 有限元法有限元法是一种将微分方程转化为代数方程组进行求解的方法。
它通过将求解域划分为有限个小区域,然后在每个小区域上选择一个试验函数来代表解,在满足边界条件的情况下,通过最小化误差的方法得到近似解。
有限元法的基本思想是将求解域划分为若干个小单元,然后在每个小单元上选择一个适当的试验函数,通过建立弱形式和加权残差方法得到代数方程组,最终通过迭代求解代数方程组得到边界条件下的解。
4.2 直接方法直接方法是通过对微分方程进行直接求解的方法,其中最常用的方法是变分法。
常微分方程的数值解法
![常微分方程的数值解法](https://img.taocdn.com/s3/m/102bbb1bc281e53a5802ff84.png)
常微分方程的数值解法在自然科学的许多领域中,都会遇到常微分方程的求解问题。
然而,我们知道,只有少数十分简单的微分方程能够用初等方法求得它们的解,多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法,逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式,通常称为近似解析方法。
还有一类近似方法称为数值方法,它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解,其中f (x, y )为x, y 的已知函数,y 0为给定的初始值,将上述问题的精确解记为y (x )。
数值方法的基本思想是:在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1,i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等,但一般取成相等的,这时na b h -=。
在这些节点上采用离散化方法,(通常用数值积分、微分。
泰勒展开等)将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来,不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1.对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时,我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从(9.2)式由于y (x 0) = y 0已给定,因而可以算出),()('000y x f x y =设x 1 = h 充分小,则近似地有:),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
常微分方程的数值解
![常微分方程的数值解](https://img.taocdn.com/s3/m/13e1bf2d10661ed9ad51f3ee.png)
f ( x, y1 ) f ( x, y2 ) L y1 y2
(其中 L 为 Lipschitz 常数)则初值问题( 1 )存 在唯一的连续解。
求问题(1)的数值解,就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便,可取 xn=x0+nh,(n=0,1,2,…), h称为步长。
(1),(2)式称为初值问题,(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组:
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法,对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法,又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式: x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程:
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn
【问题】二阶线性常微分方程边值问题的数值解法
![【问题】二阶线性常微分方程边值问题的数值解法](https://img.taocdn.com/s3/m/4f147cce5727a5e9846a614c.png)
【关键字】问题重庆三峡学院毕业论文论文题目:二阶线性常微分方程边值问题的数值解法专业:数学与应用数学年级:2004级学号:0203作者:指导老师:查中伟(教授)完成时间:2008年5月目录二阶线性常微分方程边值问题的数值解法摘要对于二阶线性常微分方程定解问题,大多数是不存在解析解的,有的方程既使存在解析解,然而解出其解也是很难办到的.尤其是工程计算中更需要的是数值解.因此,本文给出两种二阶线性常微分方程边值问题的数值解法.首先给出了利用差分法求其数值解,在介绍此方法的过程中,第一步构造了二阶线性常微分方程边值问题的差分格式(差分方程组),然后论证了该边值问题的收敛性,最后利用二阶微分的四阶紧致差分公式对第一步构造的差分格式进行精度上的改进,得到了较好的结果.接着介绍了二阶线性常微分方程边值问题的第二种数值解法——Taylor展开解法,该方法主要是先将边值问题转化为Fredholm积分方程,再经过数学处理即可得到关于近似解、近似解的一阶导数和近似解的二阶导数的线性方程组,最后利用Crammer法则解出了该二阶线性常微分方程边值问题的数值解.并且利用工程数学软件MA TLAB,给出了计算机程序,使前面两种算法在计算机上得以实现.最后给出了具体实例,分别运用以上两种解法进行求解,对这两种方法的计算精度进行了对比分析.关键词数值解;差分格式;解的收敛性;MATLABNumerical Solution for Boundary Value Problems ofthe Second-order Linear ordinary Differential EquationsKai-min Cheng(Grade 2004, mathematics and applied mathematics, School of Mathematics and Computer Science, , ChongqingWanzhou 404000)Abstract: There is no exact solution for the majority of second-order linear ordinary di fferential equations’ solution problem. Some of them even exists the exact solution, but to solve its solution is a very difficult job. Especially, we need the numerical solution urgently in Engineering Mathematics. Based on this, this paper gives two kinds of numerical solution for Boundary Value Problems of the second-order linear ordinary differential equations. Firstly, this paper gives the difference method to solve its numerical solution. In the process of introducing this method, we construct differential format of the Boundary Value Problems in the first step. Then we demonstrate the convergence of the Boundary Value Problem. Finally, we improve the accuracy for the difference format constructed in first step by the four-order differential format Secondly, this paper gives method of expansion to solve its numerical solution. In this method, we first transform the boundary value problems into Fredholm integral equation, and then can get a group of linear equations with unknowns to the approximate solution, the first order derivative of the approximate solution and the second derivative of the approximate solution after mathematical treatment, and can solve it in Crammer rule. Thirdly, we write an algorithm program by using engineering mathematical software and make above two methods realized on the computer. Finally, this paper gives a specific example and solves it with above two methods separately. This paper also analyzes the feasibility for their accuracy.Keywords: Numerical solutions; Differential format; The convergence of solutions; Matlab language0 引言当前,常微分方程的定解问题已经有很多重要结果,如解的存在性定理.在很多典型的常微分方程的解法上也有较大突破,同时也涌现出了一批较为经典的解法,如降阶法、积分变换法、变易常数法、特征方程法等方法.尽管如此,在数学领域中还存在着迄今为止还难以解出其解析解的微分方程,这就使得微分方程领域必然会产生一个新的微分方程分支——微分方程数值解法.对于较为简单的常微分方程,只需利用经典解法即可解出其解析解,边值问题也是如此.往往在实际工程中抽象出来的微分方程,其边值问题是相当复杂的,所以用求其解析解的方法来计算微分方程的边值问题往往是不适宜的,有时甚至是很难办到的.实际上,对于解微分方程,我们所要获取的或感兴趣的,往往只是一个或几个特定点上的数据,并且既使有的方程存在解析解,也并不意味着其一定能够表达成初等函数形式,如多项式函数、对数函数、指数函数、三角函数、反三角函数以及它们的有限组合形式.在利用数值计算方法理论的基础上,再辅之计算机若能很好的解决现实工程数学中的许多常微分方程的边值问题,这是非常好的.本文主要论述二阶线性常微分方程边值问题数值解的两种解法:差分法和Taylor 展开法.考虑二阶线性常微分方程的边值问题()()(),(),()y p x y q x y f x a x by a y b αβ'''++=≤≤⎧⎨==⎩(1) 1 差分法1.1 差分格式的构造对于边值问题(1)将区间[,]a b 进等距划分,分点为: 其中 b ah n-= 称0x a =与n x b =为边界点,称121,,,n x x x -为内部节点.在每个节点i x 上将,y y '''用差商近似表示.这里要求有相同的截断误差,以保证精度协调.对于内部节点,二阶导数用二阶中心差商表示,得21122()()i i i i y y y y x o h h+--+''=+,1,2,,1i n =- (2)一阶导数用一阶中心差商表示,得 211()()2i i i y y y x o h h+--'=+,1,2,,1i n =- (3)假设()i i y y x =,则1122()i i i i y y y y x h +--+''≈;11()2i i iy y y x h+--'≈ 于是方程(1)的差分方程:1111222i i i i i i i i i y y y y y p q y f h h+-+--+-++=,1,2,,1i n =- (4)其中()i i p p x =,()i i q q x =,()i i f f x =.称(4)为(1)的中心差分格式.1.2 差分格式(4)的收敛性作适当变换可以消除线性方程(1)中的一阶导数项.事实上,可令()2()xap t dt x e μ-⎰=,再作()()()y x x z x μ=,将之代入(1)得所以,不妨仅就缺少一阶导数项的方程来讨论.对于边值问题:()()(),()y q x y f x y a y b αβ''+=⎧⎨==⎩ (5) 这里假定()0q x ≤,其对应的差分问题是:11202,1,2,,1,i i i i i i n y y y q y f i n h y y αβ+--+⎧+==-⎪⎨⎪==⎩(6)下面讨论差分问题(6)的可解性.由于(6)式是关于变量(0,1,2,,)i y i n =的线性方程组,要证明它的解存在唯一,只要证明对应的齐次方程组只有零解.为此,引进文献[5]的极值原理:引理1(极值原理) 对于一组不全相等的数(0,1,2,,)i y i n =,记其中 0(1,2,,1)i q i n ≤=-如果()0(1,2,,1)i l y i n ≥=-,则i y 的正的最大值只能是0y 或n y ;如果()0i l y ≤(1,2,,1)i n =-,则i y 的负的最小值只能是0y 或n y .证明 用反证法 考察()0i l y ≥的情形,设(11)m y m n ≤≤-是正的最大值,即0max 0m i i ny y M ≤≤==>,且1m y -和1m y +中至少有一个小于M ,此时有:由于0,0m q M ≥>,故由上式推出()0m l y <,此与原假设矛盾. 此外,()0m l y ≤的情形可类似地进行讨论,证毕. 利用引理1的结论有:定理1 差分问题(6)的解存在并且是唯一的. 证明 只要证明对应的其次方程组 只有零解,由于这里()0(1,2,,1)i l y i n ==-,由极值原理知,i y 的正的最大值和负的最小值只能是0y 或n y ,而按边界条件00n y y ==,故所有()0(0,1,2,,)i l y i n ==全为零,证毕.下面运用极值原理论证差分方法的收敛性并估计误差.定理2 设i y 是差分问题(6)的解,而()i y x 是边值问题(6)的解()y x 在节点i x 处的值.则截断误差()i i i e y x y =-有估计式:2()96i M b a e h -≤(7) 式中 (4)max ()a x bM yx ≤≤=证明 由Taylor 展开式,易得2(4)111120()2()()()(),12(),()i i i i i i i i i i n y x y x y x h q y x f y x h y x y x ξξξαβ+--+⎧-++=+≤≤⎪⎨⎪==⎩(8) 将(8)与(6)相减,知误差()i i i e y x y =-满足2(4)11202()()120i i i i i i i n e e e h l e q e y h e e ξ+-⎧-+=-=⎪⎨⎪==⎩(9) 式中的i ξ一般不知道的,讨论下列差分问题211202()120i i i i i i n h l q M h εεεεεεε+-⎧-+=+=-⎪⎨⎪==⎩(9-1) 式中(4)max ()a x bM yx ≤≤=.首先证明(9)和(9-1)两个差分问题的解存在下列关系:i i e ε≤ (0,1,2,,)i n = (10)事实上,由于22(4)()()()1212i i i h h l M y l e εξ=-≤-=-,故有 又 00000;0n n n n e e e e εεεε-=-=+=+=利用引理1知0,0i i i i e e εε-≥+≥,即(10) 式成立. 我们进一步考察211202()120i i i i n h l M h ρρρρρρ+-⎧-+==-⎪⎨⎪==⎩(11) 这里()0i i i i l q ρεε-=≤,又0000n ρερε-=-=,故由引理1(注意0i q =时()i l y 就是()i l y ),知0i i ρε-≥,即i i ερ≤,于是有然而i ρ是容易求出的. 事实上,可以先求解差分问题(9)所对应的边值问题得 2()()()24h x M x a b x ρ=-- 容易验证()i i x ρρ=是差分问题(11)的解,注意到()x ρ在点2a bx +=达到最大值 因此有估计式(7),证毕.根据估计式(7)知,当0h →时有0i e →,这表明差分问题(6)是收敛的. 又因为(4)式是含有1n +个未知数(0,1,2,,)i y i n =的线性方程组,方程的个数1n -,要使方程组(4)有唯一解,还需要有两个边值条件0,n y y αβ==,它们和(4)一起构成三对角方程组.通过LU 分解,采用追赶法即可解出(4)(见3.1)1.3 差分格式(4)的改进从1.2的讨论可以得知,在1.1构造的差分格式是直接的中心差分格式,其截断误差是2()o h 即有二阶精度.基于中心差分格式的分析,直接利用二阶微分的四阶紧致差分公式,我们得到了数值求解二阶线性常微分方程边值问题的一种四阶精度的差分格式.1.3.1 改进后的差分格式的推导为了便于推导,可将方程(1)改写为()()()y f x p x y q x y '''=-- (12) 为了使格式具有更高得精度,利用[11]中的四阶差分公式2422()1/12x x y y o h h δδ⎛⎫''=+ ⎪+⎝⎭(13) 其中,2x δ为二阶中心差分算子21122i i i x i hφφφδφ+--+=(14)i φ可以表示i y 、i p 、i q 及i f 等.将(13)代入(12)得:2422()()()()()()1/12x i i i i i i x y f x p x y x q x y x o h h δδ⎛⎫'=--+ ⎪+⎝⎭(15) 即2224()()12x ii i i i i x i i i i i h y f p y q y f p y q y o h δδ''=--+--+ (16)而(2)即为22()x xx o h δφφ=+ (17)联立(16)(17)知224()()12x ii i i i i i i i i i h y f p y q y f p y q y o h δ''''=--+--+ (18) 即224(22)()12x ii i i i i i i i i i i i i i i i i i h y f p y q y f p y p y p y q y q y q y o h δ''''''''''''''''''=--+------+ (19) 显然要使(19)具有四阶精度,必须对i y '也进行四阶离散,于是利用Taylor 公式知2411()26i i i i y y h y y o h h +--''''=-+ (20) 又由(12)可得()()y y f py qy f p y py q y qy '''''''''''''''==--=---- (21) 于是有2411()()26i i i y y h y f p y py q y qy o h h +--''''''''=-----+ (22) 再将(2)、(3)、(21)及(22)代入(19),经整理并略去高阶项可得 11i i i i i i i y y y g γηλ-+++= (23)其中1111112()2312481224242412i i i i i i i i i i i i i i p p p p p p p p q q q hp q h h h γ+-+-+--+--=--+-+-+ (24) 21111112()5226246612i i i i i i i i i i i p hp q q q p p q q q h h η+-+-+----+=--++-- (25)1111112()3212481224242412i i i i i i i i i i i i i i p p p p p p p p q q q h p q h h h λ+-+-+-----=+++++-+ (26) 1122562424i i i i i i hp hp g f f f +-+-=++ (27)故差分方程(23)即为所要推导的四阶精度格式,其截断误差为4()o h .1.3.2 差分格式(23)稳定性分析现对改进后的差分格式(23)进行稳定性分析.首先引用文献[12]对差分方程的稳定性的定义. 定义1(差分算子的稳定性定义) 如果对于充分小的网格步长h ,线性差分算子h L 对任何离散函数(0,1,2,,)i v i n =均存在不依赖于h 的正常数M ,使得01max(,)max ,0i n h k k nv M v v L v i n ≤≤≤+≤≤ (28)则称差分算子h L 是稳定的.为了论证差分格式(23)的稳定性,先引入文献[13]的一个引理.引理2 假设h L 是正型差分算子,则有01111max max(,)max k n h k k nk n v v v M L v ≤≤≤≤-≤+ (29)定理3 设h L 是由(23)定义的差分算子,即11,11h i i i i i i i L y y y y i n γηλ-+≡++≤≤- (30) 其中,,,i i i γηλ由(24)—(26)确定,如果假定0i q ≤,并且有11()0i i i p p p -+-> , 11i n ≤≤- (31) 1111()min2102i i i i i i p p p h q q q -++--<+-,11i n ≤≤- (32)成立,则有(30)确定的差分算子h L 是稳定的.证明 由于,,i i i γηλ在条件(31)和(32)下满足0,0,0i i i γηλ><> (33)所以,由(30)所确定的差分算子h L 是正型的.根据引理2可以得知,存在不依赖于h 的正常数1M ,使得01111max max(,)max i n h i i ni n y y y M L y ≤≤≤≤-≤+ (34)其中,0(),()n y y a y y b αβ====.若令01max(,)M y M =,则由(34)知01max(,)max ,0i n h i i ny M y y L y i n ≤≤≤+≤≤ (35)则由差分算子稳定型的定义1得,差分算子h L 是稳定的. 证毕.定理3表明查分方程(23)的格式是稳定的,且在条件(31)和(32)满足的前提下,是对角占优()i i i ηγλ≥+的.因此,它与边值条件0(),()n y y a y y b αβ====构成的三对角方程组也可采用追赶法进行求解.2 Taylor 展式法2.1 一个Fredholm 积分方程的推导在这里可以假设()p x 、()q x 、()f x 在[,]a b 上均连续可微. 为了方便计算,现给出一类新函数的定义: 定义2 函数()()(0,1,2,)i x i γ=1()(1)(0)1()()(),1(1)!()()x i i i a x x t t dt i i x x γγγγ--⎧=-≥⎪-⎨⎪=⎩⎰ (36)其中()x γ可以为()f x 或()p x .由定义2知()00()()!()()(),0xi i i ax t p t dt i p x p x x x i '-=--≥⎰(37)现对方程(1)两边同时从a 到x 积分两次,并应用边值条件和分部积分整理得 ()(,)()()xay x V x t y t dt g x +=⎰ (38)其中(,)()()[()()]V x t p t x t q t p t '=+--;(2)()()[()()]()g x f x p a y a x a αα'=+++-方程(38)含有未知的常数()y a ',可应用(1)中的第二个边值条件()y b β=就能消除()y a ',最后导出下面的Fredholm 积分方程:()(,)()()bay x k x t y t dt h x +=⎰(39)其中11()()(,),(,)()()(,),b x b a V a t t x k x t a x b a V b t t x--⎧--≤=⎨-->⎩;1(2)(2)()()()()[()]h x f x x a b a f b ααβ-=+---+- 2.2 Taylor 展开解法设未知函数()y t 在[,]a b 上1n +阶连续可微,为了解出方程(39),我们可以将未知函数()y t 用 Taylor 公式在t x =处展开得:()1()()()()()()(,)!n n n y t y x y x t x y x t x R t x n '=+-++-+ (40) 其中(,)n R t x 为Lagrange 余项,即(1)11(,)()(),(1)!n n n R t x y t x x t n ξξ++=-≤≤+设M 为(1)()n yx +在[,]a b 上的最大值,于是得Lagrange 余项(,)n R t x 有界,即1(,)()(1)!n n MR t x b a n +≤-+ (41)特别地,若未知函数()y x 是一个次数不超过n 的多项式,则(,)0n R t x ≡.实际上,如果()y x 在[,]a b 上有任意阶的导数,只要n 充分大,余项(,)n R t x 可以忽略不计,为方便起见,本文此处仅考虑二阶的Taylor 展开式(更高阶的解法类似于二阶),即(2)21()()()()()()2y t y x y x x t y x t x '≈--+- (42) 现把 (42)代入(39)得000102[1()]()()()()()()x y x x y x x y x h x '''+∆+∆+∆= (43)其中 0(1)()()(,),0,1,2!ibi i ax x t k x t dt i i -∆=-=⎰在方程(39)两边对x 求导得 ()(,)()()bay x k x t y t dt h x ''+=⎰(44)其中 11()(,),(,)()(,),b a V a t t xk x t b a V b t t x--⎧--≤=⎨-->⎩ 现将(42)代入(44)得101112()()[1()]()()()()x y x x y x x y x h x ''''∆++∆+∆= (45)其中 1(1)()()(,),0,1,2!ibi i ax x t k x t dt i i -∆=-=⎰所以,现在方程(1)、(43)以及(45)恰好构成一个关于()y x 、()y x '和()y x '的线性方程组,应用 Crammer 法则即可解出(1)的近似解12()y x σσ=(46) 其中 010211112()()()()()()()()1h x x x h x x x f x p x σ∆∆'=∆∆,0001022101112()()()()()()()()1x x x x x x q x p x σ∆∆∆=∆∆∆由于(46)是利用Taylor 二阶展开得到的,所以(46)称为(1)的二阶近似解.如果解是二次多项式,易知二阶近似解退化成精确解,推而广之,相应的n 阶近似解对n 次多项式也是精确的.3 数值解的算法3.1 追赶法三对角线方程组的一般表示方法Ax D = (47) 其中A 为n n ⨯三对角矩阵,B 为n 阶列向量,即11222111iiin n n nn b c a b c A a b c a b c a b ---⎛⎫ ⎪ ⎪ ⎪⎪=⎪ ⎪ ⎪ ⎪ ⎪⎝⎭,121i n n x x x x x x -⎛⎫ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭,121i n n d d D d d d -⎛⎫ ⎪⎪⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎪⎝⎭并且A 满足11(1)0,0 (2),(0,2,,1)n n i i i i i b c •b c b a c a c i n ⎧>>>>⎪⎨≥+≠=-⎪⎩ (48) 设A 为满足(48)的n 阶三对角阵,则有唯一三角分解A LU =,其中L 为下三角阵,U 为单位上三角阵,即有11112222211111iiii i in nn nn b c a a b c r a A a b c r a a b r a βββ-⎛⎫⎛⎫⎛⎫⎪⎪⎪ ⎪ ⎪⎪ ⎪ ⎪⎪==⎪ ⎪⎪ ⎪ ⎪⎪ ⎪ ⎪⎪⎪ ⎪⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭(49)由矩阵乘法, 得1)11111111,,/b c c b ααββ===; 2)11,,(2,,)i i i i i i i i a r b r a i n αββ--==+==3),(2,,1)i i i c i n αβ==-于是,得到解(1)的追赶法计算公式 (1) 分解计算公式A LU =:1111//(),(2,,1)i i i i i c b c b a i n βββ-=⎧⎨=-=-⎩ (50)(2) 求解Ly D =递推公式:11111/()/(),(2,,)i i i i i i i y d b y d a y b a i n β--=⎧⎨=--=⎩ (51)(3) 求解Ux y =递推公式:1,(1,,2,1)n n i i i i x y x y x i n β+=⎧⎨=-=-⎩ (52)现将计算121n βββ-→→→及12n y y y →→→的过程称为追的过程,计算方程组解11n n x x x -→→→的过程称为赶的过程,该算法称为追赶法.定理4 对三对角线方程组(47),其中A 满足(48)式,则由追赶法计算公式得到{},{}i i αβ满足: (1) 01,(1,2,,1)i i n β<<=-(2) 0i i i i i i c b a b a α<≤-<<+,(2,,1)i n =-(3) 0n n n n n b a b a α<-<<+分析 (1)要证01i β<<,只要证1||||i i i i i a b a c β-=->,而此式中含有1i β-,因此用归纳法证明.(2)、(3)只要用三角不等式即可证得.但值得注意的是,条件(48)是充分的,但条件并非必要的,三对角线上不能有零元素,条件太苛刻,于是条件可以作适当改变,使追赶法仍可行 ,如改为:证明 (1)显然 现归纳法假设101i β-<<,下证: 01i β<<事实上,11||||||||||||||||i i i i i i i i i i a b a b a b a c ββ--=-≥->-≥,则0||||1ii ic βα<=<,(1,2,,1)i n =-(2) 因为 1i i i i a b a β-=-故 11||||||||i i i i i i i i i a b a b a b a ββ--=-≤+<+ 再有 ||||||||0i i i i a b a c >-≥>,(2,,1)i n =-(3) 由于1n n n n b a αβ-=-所以 11,(01)n n n n n n n b a b a αββ--≤+≤+<<,110,(01)n n n n n n n b a b a αββ--≥->-><<.证毕由定理4说明追赶法计算公式中不会出现中间结果数量级巨大增长和相应的舍入误差的严重累积,即追赶法计算公式对于舍入误差是稳定的.3.1.2 差分格式(23)的算法步骤1) 由1的讨论可以得知,由差分格式(23)以及边界条件0(),()n y y a y y b αβ====构成的方程组Ay g = (53) 是满足追赶法的条件(48)的. 其中,这里1111111001i iin n n A γηλγηλγηλ---⎛⎫⎪⎪⎪ ⎪=⎪ ⎪⎪ ⎪ ⎪⎝⎭,012n y y y y y ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭,012n g g g g g ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭,,,,i i i i g γηλ由(24)—(27)确定2) 根据(50)的递推公式 11111/0/10/()/(),(2,,1,)i i i i i i i i i c b c b a i n n βββληγβ--===⎧⎨=-=-=-⎩ (54)可以解出,(1,2,,1,)i i n n β=-3) 根据(51)的递推公式00001/()/(),(1,2,,)i i i i i i i z g g z g z i n ηγηγβ-==⎧⎨=--=⎩ (55)可以解出,(0,1,2,,1,)i z i n n =-4) 根据(52)的递推公式11,(1,,2,1,0)n ni i i i y z y z y i n β++=⎧⎨=-=-⎩ (56)即可解出()012,,,,Tn y y y y y = (57)故(57)即为二阶线性常微分方程边值问题(1)之差分格式(23)的数值解.3.2 Taylor 展开法的算法步骤由(46)式知,要求()y x ,只需求出1σ与2σ,而(1,2)i i σ=完全由(),(),()f x p x q x 和(0,1;0,1,2)ij i j ∆==以及()h x 及其导数()h x '确定,因为(),(),()f x p x q x 是由(1)唯一确定的,所以,以下步骤主要围绕求(0,1;0,1,2)ij i j ∆==以及()h x 及其导数()h x '来进行.1) 确定()h x 及其导数()h x ' 由(39)知1(2)(2)()()()()[()]h x f x x a b a f b ααβ-=+---+- (58) 所以要确定()h x ,首先得要确定(2)()f x ,由()()i x γ的定义知1()(1)(0)1()()(),1(1)!()()x i i i a f x x t f t dt i i f x f x --⎧=-≥⎪-⎨⎪=⎩⎰ (59)则,根据(59)作两次迭代即可求出(2)()f x ,将之代入(58)即可求出()h x ,从而解出()h x '.2) 确定(0,1;0,1,2)ij i j ∆== 由ij ∆的定义知0(1)()()(,),0,1,2!ibii a x x t k x t dt i i -∆=-=⎰;1(1)()()(,),0,1,2!ibi i ax x t k x t dt i i -∆=-=⎰(60)其中(,)k x t 和(,)k x t 分别见(39)和(44)式.所以,根据(60)以及(39)、(44)就可以确定(0,1;0,1,2)ij i j ∆==. 3)确定近似解()y x在1)和2)的基础上,由(46)即可确定二阶线性常微分方程边值问题(1)的近似解()y x .值得说明的是,以上步骤尽管较为清晰,但计算是相当繁琐的,尤其是几个定积分,所以针对几个复杂的定积分,可以适当采用数学软件Mathematica5.0进行处理.4 两种数值解法的计算机实现4.1差分格式(23)的计算机实现步骤根据1) 利用数学软件Mathematica5.0计算由(24)—(27)确定的,,,i i i i g γηλ; 2) 在Matlab 中建立一个名为chase.m 的文件,如下: function y=chase (a,b,c,g) %定义函数chase n=length(b);printf('追赶法\n'); if n-1==length(a) for i=n-1:-1:1 a(i+1)=a(i); endend%将a设置为n维向量c(1)=c(1)/b(1); f(1)=f(1)/b(1);for i=2:n-1b(i)=b(i)-a(i)*c(i-1);c(i)=c(i)/b(i);g(i)=(g(i)-a(i)*g(i-1))/b(i);endg(n)=(g(n)-a(n)*g(n-1))/(b(n)-a(n)*c(n-1));for i=n-1:-1:1g(i)=g(i)-c(i)*g(i+1);endy=f;或function g = Chase(A,g)[n,n] = size(A);printf('追赶法\n');L = tril(A);U = triu(A,1) + eye(n,n);L(1,1) = A(1,1);for k = 1:n-1L(k+1,k) = A(k+1,k);U(k,k+1) = A(k,k+1)/L(k,k);L(k+1,k+1) = A(k+1,k+1) - A(k+1,k)*U(k,k+1); enddisp(L);disp(U);b(1) = A(1,1);for k = 2:na(k) = A(k,k-1);b(k) = A(k,k);c(k-1) = A(k-1,k);B(k-1) = U(k-1,k);endz(1) = g(1)/b(1);for k = 2:nz(k) = (g(k)-a(k)*z(k-1))/(b(k)-a(k)*B(k-1)); enddisp(z);y(n) = z(n);for k = n-1:-1:1y(k) = z(k) - B(k)*y(k+1);endf = y;3) 求解线性方程组(53)Ay g =,在Matlab 中编写以上的chase.m 的M 文件,依次输入数据如下即可:>> A=[…;…;…]; >> a=[…..]; >> b=[…..]; >> c=[…..]; >> g =[…...];>> y=chase (a,b,c,g)其中A=[……]中输入(53)中A 的元素,a=[…..]中输入(1,2,,)i i n γ=,b=[…..]中输入(0,1,2,,)i i n η=,c=[…..]中输入(0,1,2,,1)i i n λ=-,g =[…...]中输入(0,1,2,,)i g i n =.4.2 Taylor 展式法的计算机实现步骤基于3.2的算法步骤,现给出更为具体的计算机操作步骤: 1) 计算()h x 及其导数()h x '首先计算(2)()f x :因为11(2)(1)1()()()[()()]1!x x x aa a f x x t f t dt x t f t dt =-=-⎰⎰⎰,所以1(2)(1)1()()()[()()]1!x x t aa a f x x t f t dt x t f u du dt =-=-⎰⎰⎰ (61)再直接在Mathematica5.0的编辑栏中输入[()()]xtaax t f u du dt -⎰⎰即可算出(2)()f x ;然后由(58)式,直接在Mathematica5.0中输入1(2)(2)()()()[()]f x x a b a f b ααβ-+---+-,即为()h x ;最后在Mathematica5.0输入D[()h x ,x ],即可得到()h x '.2) 计算(0,1;0,1,2)ij i j ∆== 由(60)式()(()()(()()))]bxa x p tb t q t p t dt '+-+--⎰ (62)同理 011()[()()(()()(()()))()x ax x t b x p t a t q t p t dt b a '∆=---+---⎰()()(()()(()()))]bxx t a x p t b t q t p t dt '+--+--⎰(63)2()()(()()(()()))]bxx t a x p t b t q t p t dt '+--+--⎰(64)(()()(()()))]bxp t b t q t p t dt '++--⎰(65)()(()()(()()))]bxx t p t b t q t p t dt '+-+--⎰ (66)2()(()()(()()))]bxx t p t b t q t p t dt '+-+--⎰ (67)所以,直接分别将(62)—(67)输入Mathematica5.0中运行后即可得到(0,1;0,1,2)ij i j ∆==3)计算近似解()y x先计算12,σσ,在Mathematica5.0中分别输入01021112()()()()()()()()1h x x x Det h x x x f x p x ∆∆⎛⎫ ⎪'∆∆ ⎪ ⎪⎝⎭和000102101112()()()()()()()()1x x x Det x x x q x p x ∆∆∆⎛⎫ ⎪∆∆∆ ⎪ ⎪⎝⎭,即可算出12,σσ.再由(46)式,可以算出(1)的近似解()y x .5 实例分析例 考虑下列边值问题2222(1)[1(1)]cos()(1)sin()1cos()(0)2,(1)1x x x y x y e y x e x x x e x y y e πππππ--'''⎧-+-=-+-++--⎨==-⎩(68) 显然,问题(68)的解析解为cos()x y e x π=+.下面分别对使用差分法和Taylor 展开法解边值问题(68)进行误差分析.5.1差分法的误差分析在(68)式中,分别取步长12341111,,,10203040h h h h ====,则由4.1的算法步骤,可以计算出改进后的差分格式(23)的最大误差max ()i i iE y y x =-以及收敛阶1212ln(/)/ln(/)r E E h h =. 结果见表1:Tab.1 Maximum error for problem (68) in differential format (23)改进后的差分格式(23)1/101291/20 8.04 4.00 1/30 1.59 4.00 1/40 0.504.005.2 Taylor 展式法的误差分析根据算法4.2可以解出(1)的近似解()y x (由于其表达式过于冗长,这里不便列出),此处只对所求结果与(1)的解析解在对应点处值进行对比分析,并且分析最大误差max ()i i iE y y x =-.同样对区间[0,1]进行40等分,只在分点处考虑问题,结果见表2:Tab.2 Maximum error for problem (68) in Taylor expansion method0 2 2 0 21/40 1.6119998 1.6119885 0.0000113 1/40 2.0222325 2.0222328 0.0000003 22/40 1.5768186 1.5768256 0.0000070 2/40 2.0389594 2.0389599 0.0000005 23/40 1.5436852 1.5437025 0.0000127 3/40 2.0202541 2.0162212 0.0000006 24/40 1.5131018 1.5132035 0.0000017 4/40 2.0562274 2.0562212 0.0000062 25/40 1.4855625 1.4856021 0.0000604 5/40 2.0562274 2.0562212 0.0000062 26/40 1.4615503 1.4615569 0.0000066 6/40 2.0528408 2.0528309 0.0000091 27/40 1.4415344 1.4418695 0.0000151 7/40 2.0438864 2.0438526 0.0000338 28/40 1.4259675 1.4250231 0.000044413/40 1.9065292 1.9065592 0.0000300 34/40 1.4486403 1.4485693 0.0000701 14/40 1.8730580 1.8730268 0.0000312 35/40 1.4749958 1.4749894 0.0000064 15/40 1.8376748 1.8375990 0.0000758 36/40 1.5085466 1.5085362 0.0000005 16/40 1.8008417 1.8008235 0.0000182 37/40 1.5494983 1.5494976 0.0000007 17/40 1.7630358 1.7630986 0.0000628 38/40 1.5980213 1.5980242 0.0000029 18/40 1.7247467 1.7247854 0.0000387 39/40 1.6542499 1.6542491 0.000000819/40 1.6864733 1.6865231 0.0000798 1 1.71828181.71828180 20/401.64872131.64869850.00007720.00042525.3 5.1与5.2结果的对比分析由表1容易得知,当设定步长1/40h =,改进后的差分格式在问题(62)的最大误差为60.5010-⨯,其收敛阶为4.00,说明该方法应用于问题(1)是相当适宜的.由表2数据对比分析得知,用Taylor 展开法来解边值问题(62),所得结果还是比较理想的.所得近似解在各节点处的值与解析解在相应节点处的值比较吻合.当然,如果(42)为更高阶的Taylor 展式,所得结果将更加精确,这说明用该法来解决二阶线性常微分方程边值问题是较为可行的.但通过对表1和表2综合对比不难发现,用改进后的差分格式(23)更能得到较为精确的数值解.致谢本论文是在指导老师查中伟教授的悉心指导下完成的.指导老师渊博的专业知识,严谨的治学态度,精益求精的工作作风,诲人不倦的高尚师德,严以律己、宽以待人的崇高风范,朴实无华、平易近人的人格魅力对我影响深远.不仅使我树立了远大的学术目标、掌握了基本的研究方法,还使我明白了许多待人接物与为人处世的道理.本论文从选题到完成,每一步都是在指导老师的指导下完成的,倾注了指导老师大量的心血.在此,谨向指导老师表示崇高的敬意和衷心的感谢!本论文的顺利完成,还离不开各位老师、同学和朋友的关心和帮助.在此感谢杨贤仆副教授、杜祥林教授、刘学飞副教授、冯天祥副教授等老师的指导和帮助;在四年的学习期间,得到了同班同学师兄师姐的大力帮助,在此表示深深的感谢.没有他们的帮助和支持是没有办法完成我的学士学位论文的,同窗、好友以及师生之间的友谊永远长存.参考文献[1] 查中伟. 数学物理偏微分方程[M]. 成都:西南交通大学出版社,2005[2] 李龙华等. 微分方程数值解法[M]. 北京:人民教育出版社,1980[3] 李群. 微分方程数值解法基础[M].北京:科学出版社,2003[4] 李大侃. 常微分方程数值解法[M].杭州:浙江大学出版社,1994[5] 李庆杨. 数值分析[M].武汉:华中理工大学出版社,2000[6] 冯天祥. 数值计算方法理论与实践研究[M], 成都:西南交通大学出版社,2005[7] 钟万勰. 结构动力学方程的精细时程积分法[J].大连理工大学学报1994,34(2):131-136[8] 魏毅强. 数值计算方法[M]. 北京:科学出版社,2004[8] 陆平. 线性二价常微分方程边值问题的讨论[J].华北工学院学报,2003,24(5):397-398[9] 赵秋林,戈新生. 求解二点边值问题打靶法的一种改进方法[J].上海力学,1999,20(4):451-458[10] 田振夫. 两点边值问题的一种高精度方法[J].贵州大学学报.1997,14(2):19-23[11] 傅得熏,马延文. 计算流体力学[M].北京:高等教育出版社.2002[12] Keller HB. Numerical Methods for Boundary Value problem [M], Waltham, MA; Blaisdell, 1968[13] 吴启光. 修正 Dennis格式的渐进解[J].数值计算与计算机的应用,1991,12(2):90-94[14] 刘会明.两点边值问题的一种高精度差分方法[J]. 上海理工大学学报,2005,(27):68-70[15] SN.Ha, CR.Lee. Numerical study for two-point boundary value problems using Green’s function Compute, Math, Apple 2002, 44:1599-1608[16] 钟献词,李显方. 第二类Fredholm积分方程的泰勒解法[J].数学理论与应用,2004(24),21-24[18] 楼顺天等. MATLAB程序语言设计[M].西安:西安电子科技大学,2003[19] 何仁斌. MATLAB工程计算及应用[M].重庆:重庆大学出版社,2001[20] 张智星. MATLAB程序语言设计与应用[M].北京:清华大学出版社,2002[21] 裘宗燕. Mathematica数学软件系统的应用及其程序设计[M].北京:北京大学出版社,1994此文档是由网络收集并进行重新排版整理.word可编辑版本!。
常微分方程边值问题的数值解法
![常微分方程边值问题的数值解法](https://img.taocdn.com/s3/m/a8eda96f51e79b8969022694.png)
……
……
……
3.6896236 3.6865656 -3.058*10-3
4.5635316 4.5612288 -2.303*10-3
5.5854269 5.5841425 -1.284*10-3
6.7725887 6.7725887 0
数值计算方法
对区间[a,b]作等距分划: x j a jh( j 0,1,2,...n)
h b a。由数值微分公式 n
y(x j )
y(x j1) y(x j1) 2h
h2 6
y(3) ( j )
y(x j )
y(x j1) 2 y(x j ) h2
y(x j1)
h2 12
y(4) ( j )
其中, j , j (x j1, x j1)。
差分方程的建立
代入y p(x) y q(x) y r(x), x [a,b]得差分方程:
y j1 2 y j h2
y j1
pj
y j1 y j1 2h
qj yj
rj ( j 1,2,...n 1)
这是求y j ( j 0,1,2,...,n)的n 1个方程,还缺的两个方
程由边界条件补出。
差分方程的建立
对于第一类边界条件:y0 , yn ,即已给出
两个未知量的解, 这时整理后有
b1 c1 a2 b2 c2 y1 d1 a1 源自y2d2an2
bn2
cn2
yn2
dn2
an1 bn1 yn1 dn1 cn1
其中a
解 : 此方程的解析解为y x2 x3 1 x4 x2 ln x。 2
例题
现取步长h 0.1,此时p(x) 2 , x
常微分方程边值问题的解法
![常微分方程边值问题的解法](https://img.taocdn.com/s3/m/4cf9897ab207e87101f69e3143323968001cf446.png)
常微分方程边值问题的解法常微分方程是描述自然科学、工程技术和经济管理等领域中各种变化规律的一个基础理论。
而边值问题是求解一些微分方程的重要问题之一,涉及到数学、物理、化学等多个领域。
在本文中,我们将讨论常微分方程边值问题的解法。
1. 边值问题的定义在微分方程解的过程中,边值问题(Boundary Value Problem, BVP)是指在区间 $[a,b]$ 上求解微分方程的解,同时已知$y(a)=\alpha$,$y(b)=\beta$ 的问题。
边值问题是对初值问题(Initial Value Problem, IVP)的一种自然延伸,在一定范围内对变量的取值进行限制,使得解的可行域更为明确。
举例来说,对于经典的二阶线性微分方程$$ y''+p(x)y'+q(x)y=f(x), \quad a<x<b $$ 如果边界条件是$y(a)=\alpha$,$y(b)=\beta$,则这个微分方程就是一个边值问题。
2. 常用解法对于一般的常微分方程边值问题,没有通用的方法可以求出其解析解,必须采用一些数值计算的方法进行求解。
常用的边值问题的解法大致有以下几种:(1)求解特殊解的方法这种方法常用于求解具有周期性边界条件的问题。
如果问题中的边界条件满足:$y(a)=y(b)=0$,则可以将问题转化为一个周期问题,即 $y(a+k)=y(b+k)$,其中 $k=b-a$。
这时,边值问题就变成了求解这个方程的周期解,例如,可以使用Fourier 级数来求解。
(2)变分法变分法是一种基于求解最小值的方法,可以用来求解一类线性边值问题。
其基本思路是将原问题转化为求一个积分的最小值。
对于一般的边值问题 $y''+f(x)y=g(x)$,可以构造一个变分问题:$$ \delta\int_a^b \left(y'^2-f(x)y^2-2gy\right) \mathrm{d}x=0 $$ 这个问题的解可以通过对变分问题的欧拉方程求解而得到。
打靶法
![打靶法](https://img.taocdn.com/s3/m/f55e584b2e3f5727a5e9626a.png)
用某种离散化数值步骤求出常微分方程边值问题在离散点上的近似解的方法。
各种实际问题导出不同类型的边值问题。
较简单的有二阶常微分方程两点边值问题:求函数y=y(x),x∈【α,b】,使它满足微分方程和边值条件式中ƒ、g1、g2为已知函数;α与b为两个给定的端点。
较一般地有一阶常微分方程组两点边值问题:求N个函数使其满足微分方程组和边值条件式中诸ƒn、g i是已知函数;r为给定的自然数。
有些问题因求解区间是无穷区间而被称作奇异边值问题,相应的边界条件变为对解在无穷远处渐近行为的限制,例如,要求y(x)在区间【0,)上平方可积或要求当x趋于无穷时,y(x)趋于某极限值。
还有些实际问题因要求解满足多个点上的条件而被称作多点边值问题。
近年来,对反映边界层现象的奇异摄动边值问题提出了一些新的数值解法。
此外,关于存在多个解的分歧现象数值解问题也引起人们的注意。
打靶法主要思路是:适当选择和调整初值条件,(选什么)求解一系列初值问题,使之逼近给定的边界条件。
如果将描述的曲线视作弹道,那么求解过程即不断调整试射条件使之达到预定的靶子,所以称作打靶法或试射法,此类方法的关键是设计选取初值的步骤。
对非线性边值问题可通过下列步骤求数值解:①计算初值问题的数值解y1。
若g(y1(b),y姈(b))=B,近似地满足,则y1即为所求;否则进行②。
②计算初值问题的数值解y2,若g(y2(b),y娦(b))=B近似地满足,则y2即为所求;否则令m=3进行③。
③将g(y(b),y┡(b))视为y(α)的函数,用线性逆插值法调整初值,即计算然后进行④。
④计算初值问题的数值解y m并进行判定:若b点边值条件近似地满足,则y m即为所求;否则令m+1崊m转向③继续计算直到满意为止。
特别地,若微分方程是线性的,则打靶法变成线性组合法,即根据常微分方程理论适当选取初值可得到一组线性独立解,利用它们的线性组合导出边值问题的解。
例如线性方程边值问题的数值解可通过两个初值问题数值解来实现。
打靶法
![打靶法](https://img.taocdn.com/s3/m/985eb94403d8ce2f01662318.png)
打靶法常微分方程边值问题数值解法- 正文用某种离散化数值步骤求出常微分方程边值问题在离散点上的近似解的方法。
各种实际问题导出不同类型的边值问题。
较简单的有二阶常微分方程两点边值问题:求函数y=y(x),x∈【α,b】,使它满足微分方程和边值条件式中ƒ、g1、g2为已知函数;α与b为两个给定的端点。
较一般地有一阶常微分方程组两点边值问题:求N个函数使其满足微分方程组和边值条件式中诸ƒn、g i是已知函数;r为给定的自然数。
有些问题因求解区间是无穷区间而被称作奇异边值问题,相应的边界条件变为对解在无穷远处渐近行为的限制,例如,要求y(x)在区间【0,)上平方可积或要求当x趋于无穷时,y(x)趋于某极限值。
还有些实际问题因要求解满足多个点上的条件而被称作多点边值问题。
近年来,对反映边界层现象的奇异摄动边值问题提出了一些新的数值解法。
此外,关于存在多个解的分歧现象数值解问题也引起人们的注意。
打靶法主要思路是:适当选择和调整初值条件,(选什么)求解一系列初值问题,使之逼近给定的边界条件。
如果将描述的曲线视作弹道,那么求解过程即不断调整试射条件使之达到预定的靶子,所以称作打靶法或试射法,此类方法的关键是设计选取初值的步骤。
对非线性边值问题可通过下列步骤求数值解:①计算初值问题的数值解y1。
若g(y1(b),y姈(b))=B,近似地满足,则y1即为所求;否则进行②。
②计算初值问题的数值解y2,若g(y2(b),y娦(b))=B近似地满足,则y2即为所求;否则令m=3进行③。
③将g(y(b),y┡(b))视为y(α)的函数,用线性逆插值法调整初值,即计算然后进行④。
④计算初值问题的数值解y m并进行判定:若b点边值条件近似地满足,则y m即为所求;否则令m+1崊m转向③继续计算直到满意为止。
特别地,若微分方程是线性的,则打靶法变成线性组合法,即根据常微分方程理论适当选取初值可得到一组线性独立解,利用它们的线性组合导出边值问题的解。
常微分方程边值问题
![常微分方程边值问题](https://img.taocdn.com/s3/m/a1687787ec3a87c24028c462.png)
常微分方程边值问题求常微分方程满足给定边界条件的解的问题。
亦即,设常微分方程为:对区间I上的点α1,α2,…,αk及值y(αi),y┡(αi),…,y(n-1)(αi)(i=1,2,…,k,k>1),给定了一些条件,求此方程在I上的满足这些条件的解的问题。
这些条件称为边界条件,诸αi及y(αi)、y┡(αi)、…、y(n-1)(αi) 称为边值或边界值。
当k=2,α1、α2是区间I的端点时,称为两点边值问题。
边值问题的提出和发展,与流体力学、材料力学、波动力学以及核物理学等密切相关;并且在现代控制理论等学科中有重要应用。
因为常微分方程可以解析求解的类型甚少,所以求边值问题的解也是困难的。
为了适应实际问题的需求,不得不采用近似解法,这样,首先需要回答:边值问题的解是否存在?是否惟一?这就是边值问题的基本论题。
在有限区间上的边值问题两点边值问题以二阶常微分方程为例。
求二阶常微分方程(1)满足边界条件的解。
式中α、b为区间的端点,ƒ:【α,b】×R2→R是连续函数,R=(-,);αs,α▂,βs,β▂及γs(s=1,2)是给定的常数。
特别当γs=0 (s=1,2)时,(2)称为齐次边界条件。
(2)的特例有:方程(1)与(2┡)、(1)与(2″)及(1)与(2冺),所构成的边值问题分别称为第一边值问题、第二边值问题和第三边值问题。
例如,悬链线(图1)之形状是由第一边值问题所确定。
式中ρ为悬链线密度,g 为重力加速度,T为悬链线最低点张力。
又如,一端固定的细长悬梁(图2)弯曲的倾斜角φ(s)是由第二边值问题Bφ″(s)-P cosφ (s) = 0,φ(0)=0,φ┡(l)=0,所确定。
其中B为梁的刚性系数,P为自由端的铅直负载。
关于边值问题解的存在和惟一性问题,对线性方程,在理论上是容易解决的。
考虑第一边值问题(3)与(2┡),其中p、q及r是【α,b】上的连续函数,设⑶的通解(4)式中y1、y2是(3)对应的齐次方程的基本解组;Y是(3)的特解;c1、c2是任意常数、为求边值问题(3)与(2┡)的解,只要将(4)代入(2┡)来确定c1、c2。
常微分方程边值问题解法
![常微分方程边值问题解法](https://img.taocdn.com/s3/m/f00f01b9cd22bcd126fff705cc17552706225e4a.png)
常微分方程边值问题解法
常微分方程边值问题解法:
常微分方程边值问题是指在一定区间内,给定一个微分方程的初始条件和边界条件,求解微分方程的解在这个区间内满足这些条件的问题。
常见的边值问题有两种类型:Dirichlet边界条件和Neumann边界条件。
解决常微分方程边值问题的方法有很多种,下面介绍其中两种常用的方法:
1. 有限差分法:
有限差分法是利用差分近似替代微分,将微分方程转化为一组代数方程。
首先将区间离散化,将连续的函数转化为离散的数值,然后利用中心差分、前向差分或后向差分的方法,将微分方程变为代数方程组,最后利用线性代数的方法求解这个方程组。
2. 有限元法:
有限元法是将区间划分为若干个小的子区间,将微分方程转化为一组局部的代数方程组,然后将这些方程组组合成整个问题的全局方程组。
有限元法可以适用于更加复杂的边值问题,但是需要更多的计算量和更高的数学水平。
总之,常微分方程边值问题的解法有很多种,需要根据具体情况选择不同的方法。
具有积分边界条件的常微分方程边值问题的数值解
![具有积分边界条件的常微分方程边值问题的数值解](https://img.taocdn.com/s3/m/bf47d8587ed5360cba1aa8114431b90d6d858941.png)
具有积分边界条件的常微分方程边值问题的数值解
具有积分边界条件的常微分方程边值问题是一种复杂的问题,要求通过
数值方法求解这种边值问题,可以有效地求解出其解,并获得非常好的精度。
首先,我们需要了解该问题的描述。
这类问题要求求解的是一个带有积
分边界条件的微分方程组。
其标准形式是: y' = f(x,y),从x0到xn的范围,其中低端为α(α为积分边界条件),高端点为β,β也被称为另一
个边界条件。
接下来,为了解决这个问题,我们首先要建立一个数值模型。
对该模型
来说,我们需要将[x0,xn]区间划分为若干等距离的子区间,并将该区间拆
分成N个等距离的点,即x1,x2,...,xn-1,它们能够完全表示该区间的函
数的取值,它们被称为子网格点。
接着,我们将把原微分方程组分解为若干
份子问题,并引入一个方程来表示积分的解的变化,求解该方程,以得到积
分边界条件的解。
最后,我们可以使用常用的数值方法求解这一模型。
一般来说,用有限
差分法求解最贴近于实际解决问题,其优点是计算量小,准确性好,不受精
度的限制。
另外,需要注意,积分边界条件可能将确定它们对应的解,有时
尤其是对于带有多个积分边界条件的情况,需要引入改进算法来获得更准确
的结果。
总而言之,具有积分边界条件的常微分方程边值问题可以通过使用数值
解的方法得到满足边界条件的数值解,有助于更好地理解微分方程的解和积
分边界条件的作用。
常微分方程边值问题的上下解方法
![常微分方程边值问题的上下解方法](https://img.taocdn.com/s3/m/c3740847302b3169a45177232f60ddccda38e6eb.png)
常微分方程边值问题的上下解方法常微分方程边值问题的上下解方法本文主要利用上下解方法研究了几类常微分方程的边值问题,得到了许多有意义的结论.第一章简要介绍了常微分边值问题上下解方法的一些发展现状.为方便阅读,在第二章给出了所要用到的定义、定理等基础知识.在第三章中,利用反序上下解方法讨论了二阶两点边值问题解的存在性和唯一性的充分条件.给出了求最小解和最大解的迭代序列,在满足解的唯一性的条件下给出了求解的迭代序列及误差估计式.第四章沿用了与第三章类似地方法讨论了二阶三点边值问题解的存在性和唯一性的充分条件,其中0 η 1,δ 0.第五章,利用上下解方法和Leray-Schauder度理论讨论了二阶两点边值问题三个解的存在性.通过求其Green函数将所要研究的边值问题的解的问题转化为求其相应算子的不动点.最后利用Leray-Schauder度理论证明了该算子的三个不动点的存在性.第六章讨论了二阶多点边值问题三个解的存在性.其中f满足Carathéodory条件,且通过求其Green函数将所要研究的边值问题的解的问题转化为求其相应算子的不动点.若f 满足Carathéodory条件和至多线性增长条件,则利用Lebesgue控制收敛定理证明了此算子的全连续性,最后应用Leray-Schauder度理论证明了算子的三个不动点的存在性.在最后一章应用和第五章类似地方法讨论了三阶三点边值问题的三个解的存在性,其中0 δ1,0 ηb.本文举了一些例子,将所得的研究结果应用于某些具体边值问题,给出了这一类边值问题解的存在性的判断方法,具有较强的实际应用背景. 摘要2-4ABSTRACT4-7第一章绪论7-10第二章预备知识10-142.1 绝对连续函数与微积分基本定理102.2 Leray-Schauder 度的一些性质10-112.3 边值问题中的Green 函数的求法11-14第三章反序上下解条件下的二阶两点边值问题解的存在性与唯一性14-21第四章反序上下解条件下的二阶三点边值问题解的存在性与唯一性21-30第五章二阶两点边值问题三个解的存在性30-41第六章二阶多点边值问题三个解的存在性41-52第七章三阶三点边值问题三个解的存在性52-62参考文献62-65在读期间公开发表的论文和承担科研项目及取得成果65-66致谢66。
第十一章 常微分方程边值问题的数值解法
![第十一章 常微分方程边值问题的数值解法](https://img.taocdn.com/s3/m/c8249c5af01dc281e43af00c.png)
第十一章 常微分方程边值问题的数值解法工程技术与科学实验中提出的大量问题是常微分方程边值问题.本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法.11.1 引言在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程),,(y y x f y '='', b x a <<, (11.1.1)在如下三种边界条件下的定解问题: 第一种边界条件:α=)(a y , β=)(b y (11.1.2)第二种边界条件:α=')(a y , β=')(b y (11.1.2)第三种边界条件:⎩⎨⎧=-'=-'1010)()()()(b b y b y a a y a y βα, (11.1.13) 其中0 0, ,00000>+≥≥b a b a .常微分方程边值问题有很多不同解法, 本书仅介绍打靶方法和有限差分方法.11.2 打靶法对于二阶非线性边值问题()()().,,βα==≤≤'=''b y a y b x a y y x f y ,,, (11.2.1)打靶法近似于使用初值求解的情况. 我们需要利用一个如下形式问题初值解的序列:()()v a w a w b x a w w x f w ='=≤≤'='')(,,,,,α, (11.2.2)引进参数v 以近似原边界值问题的解.选择参数k v v =,以使:()()β==∞→b y v b w k k ,lim , (11.2.3)其中),(k v x w 定义为初值问题(11.2.2)在k v v =时的解,同时()x y 定义为边值问题(11.2.1)的解.首先定义参数0v ,沿着如下初值问题解的曲线,可以求出点),(αa 对应的初始正视图()()v a w a w b x a w w x f w ='=≤≤'='')(,,,,,α. (11.2.4)如果),(0v b w 不严格收敛于β,那么我们选择1v 等值以修正近似值,直到),(0v b w 严格逼近β.为了取得合适的参数k v ,现在假定边值问题(11.2.1)有唯一解,如果),(v x w 定义为初始问题(11.2.2)的解,那么v 可由下式确定:0),(=-βv b w . (11.2.5)由于这是一个非线性方程,我们可以利用Newton 法求解.首先选择初始值0v ,然后由下式生成序列),)(()),((111-----=k k k k v b dvdw v b w v v β,此处),(),)((11--=k k v b dv dwv b dv dw , (11.2.6)同时要求求得),)((1-k v b dvdw,因为),(v b w 的表达式未知,所以求解这个有一点难度;我们只能得到这么一系列的值。
(完整版)常微分方程组(边值)
![(完整版)常微分方程组(边值)](https://img.taocdn.com/s3/m/16d0045d33d4b14e852468ed.png)
常微分方程组边值问题解法打靶法Shooting Method (shooting.m )%打靶法求常微分方程的边值问题function [x,a,b,n]=shooting(fun,x0,xn,eps)if nargin<3eps=1e-3;endx1=x0+rand;[a,b]=ode45(fun,[0,10],[0,x0]');c0=b(length(b),1);[a,b]=ode45(fun,[0,10],[0,x1]');c1=b(length(b),1);x2=x1-(c1-xn)*(x1-x0)/(c1-c0);n=1;while (norm(c1-xn)>=eps & norm(x2-x1)>=eps)x0=x1;x1=x2;[a,b]=ode45(fun,[0,10],[0,x0]');c0=b(length(b),1);[a,b]=ode45(fun,[0,10],[0,x1]');c1=b(length(b),1)x2=x1-(c1-xn)*(x1-x0)/(c1-c0);n=n+1;endx=x2;应用打靶法求解下列边值问题:()()⎪⎪⎩⎪⎪⎨⎧==-=010004822y y y dxy d解:将其转化为常微分方程组的初值问题()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧==-==t dx dy y y y dx dy y dx dy x 0011221048命令:x0=[0:0.1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'r')hold on[x,y]=ode45('odebvp',[0,10],[0,2]');plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,5]');plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,8]');plot(x,y(:,1))[x,y]=ode45('odebvp',[0,10],[0,10]');plot(x,y(:,1))函数:(odebvp.m)%边值常微分方程(组)函数function f=odebvp(x,y)f(1)=y(2);f(2)=8-y(1)/4;f=[f(1);f(2)];命令:[t,x,y,n]=shooting('odebvp',10,0,1e-3)计算结果:(eps=0.001)t=11.9524plot(x,y(:,1))x0=[0:1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); hold onplot(x0,y0,’o’)有限差分法Finite Difference Methods FDM (difference.m )同上例:4822y dx y d -=⇒482211i i i i y h y y y -=+--+ 2121842h y y h y i i i =+⎪⎪⎭⎫ ⎝⎛--+- 若划分为10个区间,则: ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛--⎪⎪⎭⎫ ⎝⎛----08880842114211421142222212212222h h h h y y y y h h h h n n M M O O O函数:(difference.m )%有限差分法求常微分方程的边值问题function [x,y]=difference(x0,xn,y0,yn,n)h=(xn-x0)/n;a=eye(n-1)*(-(2-h^2/4));for i=1:n-2a(i,i+1)=1;a(i+1,i)=1;endb=ones(n-1,1)*8*h^2;b(1)=b(1)-0;b(n-1)=b(n-1)-0;yy=a\b;x(1)=x0;y(1)=y0;for i=2:nx(i)=x0+(i-1)*h;y(i)=yy(i-1);endx(n)=xn;y(n)=yn;命令:[x,y]=difference(0,10,0,0,100);计算结果:x0=[0:0.1:10];y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解plot(x0,y0,'r')hold on[x,y]=difference(0,10,0,0,5);plot(x,y,’.’)[x,y]=difference(0,10,0,0,10);plot(x,y,’--’)[x,y]=difference(0,10,0,0,50);plot(x,y,’-.’)正交配置法Orthogonal Collocatioin Methods CM构造正交矩阵函数(collmatrix.m)%正交配置矩阵(均用矩阵法求对称性与非对称性正交配置矩阵)function [am,bm,wm,an,bn,wn]=collmatrix(a,m,fm,n,fn)x0=symm(a,m,fm); %a为形状因子;m为零点数;fm为对称的权函数(0为权函数1,非0为权函数1-x^2)for i=1:mxm(i)=x0(m+1-i);endxm(m+1)=1;for j=1:m+1for i=1:m+1qm(j,i)=xm(j)^(2*i-2);cm(j,i)=(2*i-2)*xm(j)^(2*i-3);dm(j,i)=(2*i-2)*(2*i-3+(a-1))*xm(j)^(2*i-3+(a-1)-1-(a-1));endfmm(j)=1/(2*j-2+a);endam=cm*inv(qm);bm=dm*inv(qm);wm=fmm*inv(qm);x1=unsymm(n,fn); %n为零点数;fn为非对称的权函数(0为权函数1,非0为权函数1-x) xn(1)=0;for i=2:n+1xn(i)=x1(n+2-i);endxn(n+2)=1;for j=1:n+2for i=1:n+2qn(j,i)=xn(j)^(i-1);if j==0 | i==1cn(j,i)=0;elsecn(j,i)=(i-1)*xn(j)^(i-2);endif j==0 | i==1 | i==2dn(j,i)=0;elsedn(j,i)=(i-2)*(i-1)*xn(j)^(i-3);endendfnn(j)=1/j;endan=cn*inv(qn);bn=dn*inv(qn);wn=fnn*inv(qn);%正交多项式求根(适用于对称问题)function p=symm(a,m,fm) %a为形状因子,m为配置点数,fm为权函数for i=1:mc1=1;c2=1;c3=1;c4=1;for j=0:i-1c1=c1*(-m+j);if fm==0c2=c2*(m+a/2+j);%权函数为1elsec2=c2*(m+a/2+j+1);%权函数为1-x^2endc3=c3*(a/2+j);c4=c4*(1+j);endp(m+1-i)=c1*c2/c4/c3;endp(m+1)=1;%为多项式系数向量,求出根后对对称问题还应开方才是零点p=sqrt(roots(p));%正交多项式求根(适用于非对称性问题)function p=unsymm(n,fn)if fn==0r(1)=(-1)^n*n*(n+1);%权函数为1elser(1)=(-1)^n*n*(n+2);%权函数为1-xendfor i=1:n-1if fn==0r(i+1)=(n-i)*(i+n+1)*r(i)/(i+1)/(i+1);%权函数为1elser(i+1)=(n-i)*(i+n+2)*r(i)/(i+1)/(i+1);%权函数为1-xendendfor j=1:np(n+1-j)=(-1)^(j+1)*r(j);endp(n+1)=(-1)^(n+1);p=roots(p);应用正交配置法求解以下等温球形催化剂颗粒内反应物浓度分布,其浓度分布的数学模型为:⎪⎪⎪⎩⎪⎪⎪⎨⎧=====⎪⎭⎫ ⎝⎛S C C r dr dC r R C dr dC r dr d r ,10,0361222 解:(1)标准化令R r x /=,S C C y /=代入微分方程及边界条件得:⎪⎪⎪⎩⎪⎪⎪⎨⎧=====⎪⎭⎫ ⎝⎛1,10,036122y x dx dy x y dx dy x dx d x(2)离散化03611=-∑+=j N i i ji y y B1,2,1+=N j Λ(3)转化为代数方程组(以3=N 为例)⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----000036363636432144434241343332312423222114131211y y y y B B B B B B B B B B B B B B B B 因为141==+y y N ,所以整理上式得:⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡+----=⎥⎥⎦⎤⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---3636363644342414321434241333231232221131211B B B B y y y B B B B B B B B B B B B 本例中的代数方程组为线性方程组,可采用线性方程组的求解方法;若为非线性方程组则采用相应的方法求解。
常微分方程的解法及应用
![常微分方程的解法及应用](https://img.taocdn.com/s3/m/b65702e2b04e852458fb770bf78a6529647d35ea.png)
常微分方程的解法及应用常微分方程是数学中的一个重要分支,广泛应用于各个领域,例如物理学、生物学、经济学等。
本文将介绍常微分方程的解法和应用。
一、常微分方程的解法常微分方程是描述物理现象和自然现象的重要数学工具,例如天文学、电子学、量子力学、流体力学、热力学、生物学、化学等。
常微分方程主要分为初值问题和边值问题两种。
1.初值问题初值问题是指在某个初始时刻$t_0$,系统的状态已知,求在此后的任意时间$t$内该系统的状态。
其一般形式如下:$$\frac{dy}{dt}=f(y,t), \ \ \ \ y(t_0)=y_0$$其中,$y$是未知的函数,$f$是已知的函数,$y_0$是已知的常数。
2.边值问题边值问题是指在某个区间$[a,b]$内,系统的状态已知,求满足某个条件的函数$y(t)$。
其一般形式如下:$$\frac{d^2y}{dt^2}=f(y,t), \ \ \ \ y(a)=y_A, \ \ \ \ y(b)=y_B$$其中,$y_A$和$y_B$是已知的常数。
3.解法常微分方程的解法有多种方法,下面介绍比较常用的两种方法:欧拉法和四阶龙格-库塔法。
(1)欧拉法欧拉法是常微分方程求解的一种最简单的数值方法,它的基本思想是将微分方程转化为差分方程,利用差分方程求解。
假设在时间t时,y的值为$y(t)$,而在时间$t+h$时的y的值可以用下式计算:$$y(t+h)=y(t)+h\times f(y(t),t)$$其中,$f(y,t)$是微分方程的右端函数,$h$是每次迭代的步长。
(2)四阶龙格-库塔法四阶龙格-库塔法是常微分方程求解的一种较为精确的数值方法,其基本思想是采用区间加权平均法对微分方程进行求解。
四阶龙格-库塔法是由四个步骤组成,分别为:1)计算斜率$k_1=f(y_i,t_i)$2)计算斜率$k_2=f(y_i+\frac{h}{2}k_1,t_i+\frac{h}{2})$3)计算斜率$k_3=f(y_i+\frac{h}{2}k_2,t_i+\frac{h}{2})$4)计算斜率$k_4=f(y_i+hk_3,t_i+h)$将这四个斜率加权平均后即得到四阶龙格-库塔法的解式:$$y_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$二、常微分方程的应用常微分方程广泛应用于各个领域,本节将介绍三个常微分方程的应用:自然增长模型、振动模型和物理模型。
常微分方程边值问题的数值解法编程
![常微分方程边值问题的数值解法编程](https://img.taocdn.com/s3/m/24ca2a44854769eae009581b6bd97f192279bf30.png)
一、概述常微分方程边值问题是科学与工程领域中常见的数学问题之一,对于这类问题的数值解法编程具有重要意义。
本文将探讨常微分方程边值问题的数值解法编程,首先介绍常微分方程边值问题的概念及数值解法的基本原理,然后详细阐述数值解法的编程实现过程,并附上相关代码示例和实际应用案例。
二、常微分方程边值问题概述在数学中,常微分方程边值问题是指给定一个常微分方程及其在一定边界条件下的解。
常微分方程边值问题的一般形式可以表示为:y''(x) = f(x, y, y'),a ≤ x ≤ b,y(a) = α,y(b) = β其中,y''(x)表示y关于x的二阶导数,f(x, y, y')为给定函数,a和b 为已知常数,α和β为已知常数。
边值问题的目标是求解这个微分方程在给定边界条件下的解y(x)。
三、常微分方程边值问题的数值解法原理1. 有限差分法有限差分法是常微分方程边值问题数值解法中常用的一种方法。
其基本思想是将微分方程中的导数用有限差分代替,然后将区域离散化为若干个网格点,最终将微分方程转化为代数方程组。
2. 有限元法有限元法是另一种常微分方程边值问题数值解法。
其核心思想是将求解区域分割为若干个单元,然后在每个单元上建立局部逼近解,最终将微分方程转化为一个线性方程组。
3. 微分方程的常用数值解法还有龙格-库塔法、预处理子法等。
四、常微分方程边值问题的数值解法编程实现1. 有限差分法的编程实现有限差分法的核心是将微分方程中的导数用有限差分代替,然后将求解区域离散化为网格点。
以二阶常微分方程为例,可以使用中心差分格式将导数近似为差分形式,然后建立代数方程组,最终通过迭代求解得到近似解。
2. 有限元法的编程实现有限元法的编程实现相对复杂,需要对求解区域进行网格剖分,建立单元和节点之间的关系,构建刚度矩阵和载荷向量,最终通过求解线性方程组得到近似解。
有限元法的编程实现通常需要使用专业的有限元软件或者自行编写有限元求解器。
数值计算中的常微分方程边值问题
![数值计算中的常微分方程边值问题](https://img.taocdn.com/s3/m/7fcc133ebb1aa8114431b90d6c85ec3a87c28b1a.png)
数值计算中的常微分方程边值问题常微分方程边值问题是数值计算中的重要研究领域之一,涉及到许多实际应用场景。
在本文中,我们将介绍边值问题的基本概念、求解方法以及应用实例。
一、什么是常微分方程边值问题在数学中,常微分方程可以使用初始值问题或边值问题来描述。
边值问题通常涉及到一个微分方程在一些给定条件下的解,而这些条件不同于初始值问题的初始条件。
对于一个二阶微分方程,如:y''(x) + p(x)y'(x) + q(x)y(x) = f(x),a < x < b其边值问题通常包含以下条件:y(a) = α,y(b) = β也就是说,我们需要找到一个函数 y(x),在满足微分方程和给定边界条件的情况下,使得 y(x) 满足问题的要求。
二、常微分方程边值问题的求解方法常微分方程边值问题的求解方法有很多种,其中最常见的是有限差分法和有限元法。
有限差分法是将微分方程在所给定的空间和时间区间内离散化,将连续的函数转换为离散的点和线段,通过计算差分方程的差分近似来求解微分方程边值问题。
这种方法的优点是计算简单,容易实现,在工科领域中应用广泛。
例如,当我们研究一条河流的河流动力学时,我们可以通过有限差分法来模拟河流的水流和流速。
有限元法是另一种流行的求解常微分方程边值问题的方法。
有限元法将微分方程的解转换为一个包含许多小单元的有限元模型。
每个有限元都可以理解为一个简单的子部件,有限元模型通过模拟这些子部件之间的相互作用来计算微分方程的解。
有限元法的优点是可以处理非线性方程,具有较高的计算精度,例如,在工程领域中,有限元法被广泛应用于机械结构力学、热传导等问题。
三、常微分方程边值问题的应用实例常微分方程边值问题可以用来解决许多实际问题,下面我们将谈谈其中的几个应用。
1. 车辆悬架设计常微分方程边值问题可以用于汽车悬架系统的设计。
当车辆行驶在不平路面上时,悬架系统需要运作以使车辆保持平衡和稳定性。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
第十一章 常微分方程边值问题的数值解法工程技术与科学实验中提出的大量问题是常微分方程边值问题.本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法.11.1 引言在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程),,(y y x f y '='', b x a <<, (11.1.1)在如下三种边界条件下的定解问题: 第一种边界条件:α=)(a y , β=)(b y (11.1.2)第二种边界条件:α=')(a y , β=')(b y (11.1.2)第三种边界条件:⎩⎨⎧=-'=-'1010)()()()(b b y b y a a y a y βα, (11.1.13) 其中0 0, ,00000>+≥≥b a b a .常微分方程边值问题有很多不同解法, 本书仅介绍打靶方法和有限差分方法.11.2 打靶法对于二阶非线性边值问题()()().,,βα==≤≤'=''b y a y b x a y y x f y ,,, (11.2.1)打靶法近似于使用初值求解的情况. 我们需要利用一个如下形式问题初值解的序列:()()v a w a w b x a w w x f w ='=≤≤'='')(,,,,,α, (11.2.2)引进参数v 以近似原边界值问题的解.选择参数k v v =,以使:()()β==∞→b y v b w k k ,lim , (11.2.3)其中),(k v x w 定义为初值问题(11.2.2)在k v v =时的解,同时()x y 定义为边值问题(11.2.1)的解.首先定义参数0v ,沿着如下初值问题解的曲线,可以求出点),(αa 对应的初始正视图()()v a w a w b x a w w x f w ='=≤≤'='')(,,,,,α. (11.2.4)如果),(0v b w 不严格收敛于β,那么我们选择1v 等值以修正近似值,直到),(0v b w 严格逼近β.为了取得合适的参数k v ,现在假定边值问题(11.2.1)有唯一解,如果),(v x w 定义为初始问题(11.2.2)的解,那么v 可由下式确定:0),(=-βv b w . (11.2.5)由于这是一个非线性方程,我们可以利用Newton 法求解.首先选择初始值0v ,然后由下式生成序列),)(()),((111-----=k k k k v b dvdw v b w v v β,此处),(),)((11--=k k v b dv dwv b dv dw , (11.2.6)同时要求求得),)((1-k v b dvdw,因为),(v b w 的表达式未知,所以求解这个有一点难度;我们只能得到这么一系列的值。
,,,),(),(),(),(1210-⋯⋯k v b w v b w v b w v b w假如我们如下改写初值问题(11.2.2),使其强调解对x 和v 的依赖性()()v v a w v a w b x a v x w v x w x f w ='=≤≤'=''),(,),(),,(,,,,α,(11.2.7)保留初始记号以显式与x 的微分相关.既然要求当k v v =时),)((v b dvdw的值,那么我们需要求出表达式(11.2.7)关于v 的偏导数.过程如下:)),(),,(,(),(v x w v x w x v fv x v w '∂∂=∂''∂),()),(),,(,()),(),,(,(v x v w v x w v x w x w f v x v x w v x w x x f ∂∂'∂∂+∂∂'∂∂=),()),(),,(,(v x v w v x w v x w x w f ∂'∂''∂∂+又因为x 跟v 相互独立,所以当b x a ≤≤上式如下;),()),(),,(,(),()),(),,(,(),(v x vw v x w v x w x w f v x v w v x w v x w x w f v x v w ∂'∂''∂∂+∂∂'∂∂=∂''∂(11.2.8)初始条件为:1),(0),(=∂'∂=∂∂v a v w v a v w , 如果简单地用),(v x z 定义),)((v x vw∂∂,并且假定x 和v 的微分阶翻转,(11.2.8)转化为初值问题:1)(0)(),,(),,(....==≤≤''∂∂+'∂∂=a z a z b x a z w w x w fz w w x w f z ,,,,(11.2.9) 因此,对于(11.2.2)和(11.2.9)式的每次迭代需要求解两个初值问题.那么从(11.2.6)式可得:),(),(111-----=k k k k v b z v b w v v β, (11.2.10)事实上,这些初值问题很难精确求解,而将这些解近似为一个初值问题的解.同样,我们可以按以上步骤考虑对于三阶非线性边值问题的打靶法算法. 对于三阶非线性边值问题()()().)(,,,,γβα='='=≤≤'''='''b y a y a y b x a y y y x f y ,,,(11.2.11)转变形式:()()v a w a w a w b x a w w w x f w =''='=≤≤'''=''')()(,,,,,,,βα,(11.2.12)选择参数k v v =,以使:()()β=='∞→b y v b w k k ,lim , (11.2.13)其中),(k v x w '定义为初值问题(11.2.12)在k v v =时的解,同时()x y 定义为边值问题(11.2.11)的解.定义参数0v ,沿着如下初值问题解的曲线,可以求出点),(αa 对应的初始正视图()()v a w a w a w b x a w w w x f w =''='=≤≤'''=''')()(,,,,,,,βα.(11.2.14)如果),(0v b w '不严格收敛于β,那么我们选择1v 等值以修正近似值,直到),(0v b w '严格逼近β.为了取得合适的参数k v ,现在假定边值问题(11.2.11)有唯一解,如果),(v x w 定义为初始问题(11.2.12)的解,那么v 可由下式确定:0),(=-'βv b w . (11.2.15)利用Newton 法求解.首先选择初始值0v ,然后由下式生成序列),)(()),((111---'-'-=k k k k v b dvw d v b w v v β,此处),(),)((11--'='k k v b dv w d v b dv w d , (11.2.16)同时要求求得),)((1-'k v b dvw d ,因为),(v b w '的表达式未知,所以只能得到这么一系列的值。
,,,),(),(),(),(1210-'⋯⋯'''k v b w v b w v b w v b w改写初值问题(11.2.12),使其强调解对x 和v 的依赖性()(),(,),(,),(,),(,)(,)w f x w x v w x v w x v a x b w a v w a v w a v vαβ''''''='''≤≤===,,,, (11.2.17)要求当k v v =时),)((v b dvw d '的值,那么我们需要求出表达式(11.2.17)关于v 的偏导数.所以当b x a ≤≤上式如下;vw w f v w w f v w w f v x v w ∂''∂''∂∂+∂'∂'∂∂+∂∂∂∂=∂'''∂),( (11.2.18) 初始条件为:1),(0),(0),(=∂''∂=∂'∂=∂∂v a vw v a v w v a v w ,, 对于第三种边界条件, 注意到 01)(a a y'(a)a y -=, 取 s a y =)('得到 ⎪⎩⎪⎨⎧=-==s a , y'a a a y'a y x,y,y'f y )()()()(01, b x a <<, (11.2.19) 设它的解为)(x,s y , 代入第二个边界条件得到s βs F s b y βb,s y b y βb y ==-=-')() ,()(')()(00即s βs F =)(当1s ββ= 时, 即为所求的解. 这样,同第一种边界条件一样, 可以利用打靶方法求解.例11.1 利用打靶法求解两点边值问题2()0(0)0,(0.4) 1.8y y x y y '''⎧-+=⎨==⎩ 解:(1) 为应用打靶方法,需要假定初值(0)y '来先求解初值问题,取初始参数0 1.8120.40t -==-, 解2()(0)0,(0)2y y x y y '''⎧=-⎨'==⎩令1y u =,2y u '=,则上述初值问题变为1212222,(0)1,(0)2u u u u u x u ⎧'==⎪⎨'=-=⎪⎩, 由标准的R-K 方法,取0.2h =,可得0 2.540β=。
(2)再令11t =, 解1212222,(0)1,(0)1u u u u u x u ⎧'==⎪⎨'=-=⎪⎩, 取0.2h =。
仍由标准的R-K 方法,可得1 1.497 1.8β=<。
(3)令2121(1.8 1.497) 1.291.497 2.540t -=+-=-再解1212222,(0)1,(0) 1.29u u u u u x u ⎧'==⎪⎨'=-=⎪⎩, 取0.11h =,得2 1.7094β=。