第八章常微分方程数值解
- 格式:doc
- 大小:492.00 KB
- 文档页数:13
求常微分方程的数值解一、背景介绍常微分方程(Ordinary Differential Equation,ODE)是描述自然界中变化的数学模型。
常微分方程的解析解往往难以求得,因此需要寻找数值解来近似地描述其行为。
求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。
二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。
它基于导数的定义,将微分方程转化为差分方程,通过迭代计算得到近似解。
欧拉法的公式如下:$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中,$y_n$表示第$n$个时间步长处的函数值,$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数,$\Delta t$表示时间步长。
欧拉法具有易于实现和理解的优点,但精度较低。
2. 改进欧拉法(Heun方法)改进欧拉法又称Heun方法或两步龙格-库塔方法,是对欧拉法进行了精度上提升后得到的一种方法。
它利用两个斜率来近似函数值,并通过加权平均来计算下一个时间步长处的函数值。
改进欧拉法的公式如下:$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高,但计算量也更大。
3. 龙格-库塔法(RK4方法)龙格-库塔法是求解常微分方程中最常用的数值方法之一。
它通过计算多个斜率来近似函数值,并通过加权平均来计算下一个时间步长处的函数值。
RK4方法是龙格-库塔法中最常用的一种方法,其公式如下:$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程,可以通过以下步骤求解其数值解:1. 确定初值条件:确定$t=0$时刻的函数值$y(0)$。
常微分方程的数值解算法常微分方程的数值解算法是一种对常微分方程进行数值计算的方法,这可以帮助我们更好地理解和研究自然现象和工程问题。
在本文中,我们将介绍一些常用的数值解算法,探讨它们的优缺点和适用范围。
常微分方程(ODE)是描述自然现象和工程问题的重要数学工具。
然而,对于许多ODE解析解是无法求出的,因此我们需要通过数值方法对其进行求解。
常微分方程可以写作:y' = f(t, y)其中,y是函数,f是给定的函数,表示y随t的变化率。
这个方程可以写成初始值问题(IVP)的形式:y'(t) = f(t,y(t)),y(t0) = y0其中,y(t0)=y0是方程的初始条件。
解决IVP问题的典型方法是数值方法。
欧拉方法欧拉方法是最简单的一阶数值方法。
在欧拉方法中,我们从初始条件开始,并在t = t0到t = tn的时间内,用以下公式逐步递推求解:y n+1 = y n + hf (t n, y n)其中,f(t n,y n)是点(t n,y n)处的导数, h = tn - tn-1是时间间隔。
欧拉方法的优点是简单易懂,容易实现。
然而,它的缺点是在整个时间段上的精度不一致。
程度取决于使用的时间间隔。
改进的欧拉方法如果我们使用欧拉方法中每个时间段的中间点而不是起始点来估计下一个时间点,精度就会有所提高。
这个方法叫做改进的欧拉方法(或Heun方法)。
公式为:y n+1 = y n + h½[f(t n, y n)+f(tn+1, yn + h f (tn, yn))]这是一个二阶方法,精度比欧拉方法高,但计算量也大一些。
对于易受噪声干扰的问题,改进的欧拉方法是个很好的选择。
Runge-Kutta方法Runge-Kutta方法是ODE计算的最常用的二阶和高阶数值方法之一。
这个方法对定义域内的每个点都计算一个导数。
显式四阶Runge-Kutta方法(RK4)是最常用的Runge-Kutta方法之一,并已得到大量实践的验证。
常微分方程的数值解法常微分方程是研究变量的变化率与其当前状态之间的关系的数学分支。
它在物理、工程、经济等领域有着广泛的应用。
解常微分方程的精确解往往十分困难甚至不可得,因此数值解法在实际问题中起到了重要的作用。
本文将介绍常见的常微分方程的数值解法,并比较其优缺点。
1. 欧拉方法欧拉方法是最简单的数值解法之一。
它基于近似替代的思想,将微分方程中的导数用差商近似表示。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)根据微分方程的定义使用近似来计算下一个点的值。
欧拉方法的计算简单,但是由于误差累积,精度较低。
2. 改进欧拉方法为了提高欧拉方法的精度,改进欧拉方法应运而生。
改进欧拉方法通过使用两个点的斜率的平均值来计算下一个点的值。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)根据微分方程的定义使用近似来计算下一个点的值。
改进欧拉方法相较于欧拉方法而言,精度更高。
3. 龙格-库塔法龙格-库塔法(Runge-Kutta)是常微分方程数值解法中最常用的方法之一。
它通过迭代逼近精确解,并在每一步中计算出多个斜率的加权平均值。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)计算各阶导数的导数值。
(4)根据权重系数计算下一个点的值。
与欧拉方法和改进欧拉方法相比,龙格-库塔法的精度更高,但计算量也更大。
4. 亚当斯法亚当斯法(Adams)是一种多步法,它利用之前的解来近似下一个点的值。
具体步骤如下:(1)确定初始条件,即问题的初值。
(2)选择相应的步长h。
(3)通过隐式或显式的方式计算下一个点的值。
亚当斯法可以提高精度,并且比龙格-库塔法更加高效。
5. 多步法和多级法除了亚当斯法,还有其他的多步法和多级法可以用于解常微分方程。
多步法通过利用多个点的值来逼近解,从而提高精度。
而多级法则将步长进行分割,分别计算每个子问题的解,再进行组合得到整体解。
第8章 常微分方程边值问题的数值解法8.1 引 言推论 若线性边值问题()()()()()(),,(),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤⎧⎨==⎩ (8.1.2) 满足(1) (),()p x q x 和()f x 在[,]a b 上连续; (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。
求边值问题的近似解,有三类基本方法:(1) 差分法(difference method),也就是用差商代替微分方程及边界条件中的导数,最终化为代数方程求解;(2) 有限元法(finite element method);(3) 把边值问题转化为初值问题,然后用求初值问题的方法求解。
8.2 差分法8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法设二阶线性常微分方程的边值问题为(8.2.1)(8.2.2)()()()(),,(),(),y x q x y x f x a x b y a y b αβ''-=<<⎧⎨==⎩其中(),()q x f x 在[,]a b 上连续,且()0q x ≥.用差分法解微分方程边值问题的过程是:(i) 把求解区间[,]a b 分成若干个等距或不等距的小区间,称之为单元;(ii) 构造逼近微分方程边值问题的差分格式. 构造差分格式的方法有差分法, 积分插值法及变分插值法;本节采用差分法构造差分格式;(iii) 讨论差分解存在的唯一性、收敛性及稳定性;最后求解差分方程. 现在来建立相应于二阶线性常微分方程的边值问题(8.2.1), (8.2.2)的差分方程. ( i ) 把区间[,]I a b =N 等分,即得到区间[,]I a b =的一个网格剖分:011N N a x x x x b -=<<<<=,其中分点(0,1,,)i x a ih i N =+=,并称之为网格节点(grid nodes);步长b a Nh -=. ( ii ) 将二阶常微分方程(8.2.2)在节点i x 处离散化:在内部节点(1,2,,1)i x i N =-处用数值微分公式2(4)1112()2()()()(),12i i i i i i i i y x y x y x h y x y x x h ξξ+---+''=-<< (8.2.3)代替方程(8.2.2)中()i y x '',得112()2()()()()()()i i i i i i i y x y x y x q x y x f x R x h +--+-=+,(8.2.4) 其中2(4)()()12i i h R x y ξ=. 当h 充分小时,略去式(8.2.4)中的()i R x ,便得到方程(8.2.1)的近似方程1122i i i i i i y y y q y f h +--+-=,(8.2.5)其中(),()i i i i q q x f f x ==, 11,,i i i y y y +-分别是11(),(),()i i i y x y x y x +-的近似值, 称式(8.2.5)为差分方程(difference equation),而()i R x 称为差分方程(8.2.5)逼近方程(8.2.2)的截断误差(truncation error). 边界条件(8.7.2)写成0,.N y y αβ==(8.2.6)于是方程(8.2.5), (8.2.6)合在一起就是关于1N +个未知量01,,,N y y y ,以及1N +个方程式的线性方程组:2211212211222111(2),(2),1,2,,1,(2).i i i i i N N N N q h y y h f y q h y y h f i N y q h y h f αβ-+----⎧-++=-⎪-++==-⎨⎪-+=-⎩(8.2.7)这个方程组就称为逼近边值问题(8.2.1), (8.2.2)的差分方程组(system of difference equations)或差分格式(difference scheme),写成矩阵形式2211122222223332222222111(2)11(2)11(2)11(2)11(2)N N N N N N y q h h f y q h h f y q h h f y q h h f y q h h f αβ------⎡⎤⎡⎤-+-⎡⎤⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦. (8.2.8)用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或(8.2.8), 其解01,,,N y y y 称为边值问题(8.2.1), (8.2.2)的差分解(difference solution). 由于(8.2.5)是用二阶中心差商代替方程(8.2.1)中的二阶微商得到的,所以也称式(8.2.7)为中心差分格式(centered-difference scheme).( iii ) 讨论差分方程组(8.2.7)或(8.2.8)的解是否收敛到边值问题(8.2.1), (8.2.2)的解,估计误差.对于差分方程组(8.2.7),我们自然关心它是否有唯一解;此外,当网格无限加密,或当0h →时,差分解i y 是否收敛到微分方程的解()i y x . 为此介绍下列极值原理:定理8.2.1 (极值原理) 设01,,,N y y y 是给定的一组不全相等的数,设1122(),0,1,2,,i i i i i i i y y y l y q y q i N h +--+=-≥=.(8.2.9)(1) 若()0,1,2,,i l y i N ≥=, 则{}0Ni i y =中非负的最大值只能是0y 或N y ; (2) 若()0,1,2,,i l y i N ≤=, 则{}0Ni i y =中非正的最小值只能是0y 或N y .证 只证(1) ()0i l y ≥的情形,而(2) ()0i l y ≤的情形可类似证明. 用反证法. 记{}0max i i NM y ≤≤=,假设0M ≥, 且在121,,,N y y y -中达到. 因为i y 不全相等,所以总可以找到某个00(11)i i N ≤≤-,使0i y M =,而01i y -和01i y +中至少有一个是小于M 的. 此时0000000011222()2.i i i i i i i i y y y l y q y h M M M q M q M h +--+=--+<-=-因为00,0i q M ≥≥,所以0()0i l y <, 这与假设矛盾,故M 只能是0y 或N y . 证毕!推论 差分方程组(8.2.7)或(8.2.8)的解存在且唯一. 证明 只要证明齐次方程组11202()0,0,1,2,,1,0,0i i i i i i i N y y y l y q y q i N h y y +--+⎧=-=≥=-⎪⎨⎪==⎩ (8.2.10)只有零解就可以了. 由定理8.7.1知,上述齐次方程组的解01,,,N y y y 的非负的最大值和非正的最小值只能是0y 或N y . 而00,0N y y ==,于是0,1,2,,.i y i N == 证毕!利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果: 定理8.2.2 设i y 是差分方程组(8.2.7)的解,而()i y x 是边值问题(8.2.1), (8.2.2)的解()y x 在i x 上的值,其中0,1,,i N =. 则有224()(),96i i i M h y x y b a ε=-≤-(8.2.11)其中(4)4max ()a x bM yx ≤≤=.显然当0h →时,()0i i i y x y ε=-→. 这表明当0h →时,差分方程组(8.2.7)或(8.2.8)的解收敛到原边值问题(8.7.1), (8.7.2)的解.例8.2.1 取步长0.1h =,用差分法解边值问题43,01,(0)(1)0,y y x x y y ''-=≤≤⎧⎨==⎩并将结果与精确解()()2222()3434x xy x e e ee x --=---进行比较.解 因为110N h ==,()4,()3q x f x x ==, 由式(8.2.7)得差分格式221222112289(240.1)30.10.1,(240.1)30.1,2,3,,8,(240.1)30.10.9,i i i i y y y y y x i y y -+⎧-+⨯+=⨯⨯⎪-+⨯+=⨯=⎨⎪-+⨯=⨯⨯⎩0100y y ==, 00.1,1,2,,9i x ih i i =+==, 其结果列于表8.2.1.从表8.2.1可以看出, 差分方法的计算结果的精度还是比较高的. 若要得到更精确的数值解,可用缩小步长h 的方法来实现.8.2.2 一般二阶线性常微分方程边值问题的差分法对一般的二阶微分方程边值问题1212()()()()()(),,()(),()(),y x p x y x q x y x f x a x b y a y a y b y b αααβββ'''++=<<⎧⎪'+=⎨⎪'+=⎩ (8.2.12) 假定其解存在唯一.为求解的近似值,类似于前面的做法,( i ) 把区间[,]I a b =N 等分,即得到区间[,]I a b =的一个网格剖分:011N N a x x x x b -=<<<<=,其中分点(0,1,,)i x a ih i N =+=,步长b a Nh -=. ( ii ) 对式(8.2.12)中的二阶导数仍用数值微分公式2(4)1112()2()()()(),12i i i i i i i iy x y x y x h y x y x x h ξξ+---+''=-<<代替,而对一阶导数,为了保证略去的逼近误差为2()O h ,则用3点数值微分公式;另外为了保证内插,在2个端点所用的3点数值微分公式与内网格点所用的公式不同,即21112012000022212()()()(),,1,2,,1,263()4()()()(),,23()4()3()()(),.23i i i i i i i N N N N N N N N y x y x h y x y x x i N h y x y x y x h y x y x x h y x y x y x h y x y x x h ξξξξξξ+-----⎧-''''=-<<=-⎪⎪-+-⎪''''=+<<⎨⎪⎪-+''''=+<<⎪⎩(8.2.13) 略去误差,并用()i y x 的近似值i y 代替()i y x ,(),(),()i i i i i i p p x q q x f f x ===,便得到差分方程组1111221001221211(2)(),1,2,,1,2(34),2(43),2i i i i i i i i i N N N N p y y y y y q y f i N h hy y y y h y y y y h αααβββ-++---⎧-++-+==-⎪⎪⎪+-+-=⎨⎪⎪+-+=⎪⎩(8.2.14)其中(),(),(),1,2,,1i i i i i i q q x p p x f f x i N ====-, i y 是()i y x 的近似值. 整理得12021222211222121(23)42,(2)2(2)(2)2,1,2,,1,4(32)2.i i i i i i i N N N h y y y h hp y h q y hp y h f i N y y h y h αααααβββββ-+---+-=⎧⎪---++==-⎨⎪-++=⎩ (8.2.15)解差分方程组(8.2.15),便得边值问题(8.2.12)的差分解01,,,N y y y .特别地, 若12121,0,1,0ααββ====,则式(8.2.12)中的边界条件是第一类边值条件:(),();y a y b αβ==此时方程组(7.7.16)为221112112211221211112(2)(2)2(2),(2)2(2)(2)2,2,3,,2,(2)2(2)2(2).i i i i i i i N N N N N N h q y hp y h f hp hp y h q y hp y h f i N hp y h q y h f hp αβ-+------⎧--++=--⎪---++==-⎨⎪---=-+⎩(8.2.16) 方程组(8.2.16)是三对角方程组,用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.16),便得边值问题(8.2.12)的差分解01,,,N y y y .( iii ) 讨论差分方程组(8.2.16)的解是否收敛到微分方程的解,估计误差. 这里就不再详细介绍.例8.2.2取步长/16h π=,用差分法求下列边值问题的近似解,并将结果与精确解进行比较.精确解是1()(sin 3cos )10y x x x =-+. 解 因为(20)8N h π=-=,()1,()2,()cos p x q x f x x =-=-=, 由式(8.2.17)得差分格式()()()()()()()()()()()()()2122211222122216(2)216(1)216cos 16216(1)(0.3),216(1)2216(2)216(1)216cos 16,2,3,,6,216(1)2216(2)216cos 7i i i N N y yy y y i i y y πππππππππππππ-+--⎡⎤--⨯-++⨯-⎡⎤⎣⎦⎣⎦=--⨯-⨯-⎡⎤⎣⎦⎡⎤-⨯---⨯-++⨯-⎡⎤⎡⎤⎣⎦⎣⎦⎣⎦==⎡⎤-⨯---⨯-⎡⎤⎣⎦⎣⎦=()()16216(1)(0.1),ππ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪-+⨯-⨯-⎡⎤⎣⎦⎩080.3,0.1y y =-=-, 00.1,1,2,,9i x ih i i =+==, 其结果列于表8.2.2.8.3 有限元法有限元法(finite element method)是求解微分方程定解问题的有效方法之一,它特别适用在几何、物理上比较复杂的问题. 有限元法首先成功地应用于结构力学和固体力学,以后又应用于流体力学、物理学和其他工程科学. 为简明起见,本节以线性两点边值问题为例介绍有限元法.考虑线性两点边值问题()(8.3.1)(8.3.2)()()()()(),,(),(),Ly p x y x q x y x f x a x b y a y b αβ⎧''⎪=-+=≤≤⎨==⎪⎩其中1()0,()0,C [,]p x q x p a b >≥∈, ,C[,]q f a b ∈.此微分方程描述了长度为b a -的可变交叉截面(表示为()q x )的横梁在应力()p x 和()f x 下的偏差()y x .8.3.1 等价性定理记{}221C [,]()C [,],(),()a b y y y x a b y a y b αβ==∈==, 引进积分()22()()[()]()()2()()d baI y p x y x q x y x f x y x x '=+-⎰. (8.3.3)任取21()C [,]y y x a b =∈,就有一个积分值()I y 与之对应,因此()I y 是一个泛函(functional),即函数的函数. 因为这里是,y y '的二次函数,因此称()I y 为二次泛函.对泛函(8.3.3)有如下变分问题(variation problem):求函数21C [,]y a b ∈,使得对任意21C [,]y a b ∈, 均有()()I y I y ≥, (8.3.4)即()I y 在y 处达到极小, 并称y 为变分问题(8.3.4)的解.可以证明:定理8.3.1(等价性定理) y 是边值问题(8.3.1), (8.3.2)的解的充分必要条件是y 使泛函()I y 在21C [,]a b 上达到极小,即y 是变分问题(8.3.4)在21C [,]a b 上的解. 证 (充分性) 设21C [,]y a b ∈是变分问题()I y 的解;即y 使泛函()I y 在21C [,]a b 上达到极小,证明y 必是边值问题(8.3.1), (8.3.2)的解.设()x η是2C [,]a b 任意一个满足()()0a b ηη==的函数,则函数21()()()C [,]y x y x x a b αη=+∈,其中α为参数. 因为y 使得()I y 达到极小,所以()()I y I y αη+≥,即积分()22()()[()()]()[()()]2()[()()]baI y p x y x x q x y x x f x y x x dxαηαηαηαη''+=+++-+⎰作为α的函数,在0α=处取极小值()I y ,故d()0d I y ααηα=+=. (8.3.5) 计算上式,得()()()()()022(8.d()d d ()[()()]()[()()]2()[()()]d d 2()[()()]()2()[()()]()2()()d 2()()()()()()()()d .bab abaI y p x y x x q x y x x f x y x x x p x y x x x q x y x x x f x x x p x y x x q x y x x f x x x ααααηααηαηαηααηηαηηηηηη===+''=+++-+'''=+++-''=+-⎰⎰⎰3.6)利用分部积分法计算积分[][]()()()d ()()d ()()()()()()()d ()()()d ,bbaab ba abap x y x x x p x y x x p x y x x x p x y x x x p x y x x ηηηηη'''='''=-''=-⎰⎰⎰⎰代入式(8.3.6),得()0(8.3.7)d()2()()()()()()d 0.d b a I y p x y x q x y x f x x x ααηηα'=⎡⎤⎣⎦'+=-+-=⎰因为()x η是任意函数,所以必有()()()()()()0p x y x q x y x f x ''-+-≡. (8.3.8)否则,若在[,]a b 上某点0x 处有()00000()()()()()0p x y x q x y x f x ''-+-≠,不妨设()00000()()()()()0p x y x q x y x f x ''-+->,则由函数的连续性知,在包含0x 的某一区间00[,]a b 上有()()()()()()0p x y x q x y x f x ''-+->.作002200000,[,]\[,],()()(),.x a b a b x x a x b a x b η∈⎧⎪=⎨--≤≤⎪⎩ 显然2()C [,]x a b η∈,且()()0a b ηη==,但()()()()()()()d ba p x y x q x y x f x x x η⎡⎤''-+-⎢⎥⎣⎦⎰ ()00()()()()()()d 0b a p x y x q x y x f x x x η⎡⎤''=-+->⎢⎥⎣⎦⎰,这与式(8.3.7)矛盾. 于是式(8.3.8)成立,即变分问题(8.3.4)的解y 满足微分方程(8.3.1), 且(),()y a y b αβ==故它是边值问题(8.3.1), (8.3.2)的解.(必要性) 设()y y x =是边值问题(8.3.1), (8.3.2)的解,证明y 是变分问题(8.3.4)的解;即:y 使泛函()I y 在21C [,]a b 上达到极小.因为()y y x =满足方程(8.3.1),所以()()()()()()0p x y x q x y x f x ''-+≡.设任意21()C [,]y y x a b =∈,则函数()()()x y x y x η=-满足条件()()0a b ηη==,且2()C [,]x a b η∈. 于是()()[]()222222()()()()()[()()]()[()()]2()[()()]d ()[()]()[()]2()()d 2()()()()()()()()d ()[()]()[()]d baba baaI y I y I y I y p x y x x q x y x x f x y x x x p x y x q x y x f x y x xp x y x x q x y x x f x x x p x x q x x xηηηηηηηηη-=+-''=+++-+'-+-''=+-++⎰⎰⎰()()()22222()()()()()()d ()[()]()[()]d ()[()]()[()]d .bb ba a bap x y x q x y x f x x x p x x q x x xp x x q x x x ηηηηη⎡⎤'''=--+++⎢⎥⎣⎦'=+⎰⎰⎰⎰因为()0,()0p x q x >≥,所以当()0x η≠时,()22()[()]()[()]d 0bap x x q x x x ηη'+>⎰, 即()()0I y I y ->.只有当()0x η≡时,()()0I y I y -=. 这说明y 使泛函()I y 在21C [,]a b 上达到极小. 证毕!定理8.3.2 边值问题(8.3.1), (8.3.2)存在唯一解.证明 用反证法. 若12(),()y x y x 都是边值问题(8.3.1), (8.3.2)的解,且不相等,则由定理8.3.1知,它们都使泛函()I y 在21C [,]a b 上达到极小,因而12()()I y I y > 且 21()()I y I y >,矛盾!因此边值问题(8.3.1), (8.3.2)的解是唯一的.由边值问题解的唯一性,不难推出边值问题(8.3.1), (8.3.2)解的存在性(留给读者自行推导).8.3.2 有限元法等价性定理说明,边值问题(8.3.1), (8.3.2)的解可化为变分问题(8.3.4)的求解问题. 有限元法就是求变分问题近似解的一种有效方法. 下面给出其解题过程:第1步 对求解区间进行网格剖分01,i n a x x x x b =<<<<<=区间1[,]i i i I x x -=称为单元,长度1(1,2,,)i i i h x x i n -=-=称为步长,1max i i nh h ≤≤=. 若(1,2,,)i h h i n ==,则称上述网格剖分为均匀剖分.给定剖分后,泛函(8.3.3)可以写成()22()()[()]()()2()()d baI y p x y x q x y x f x y x x '=+-⎰()12211()[()]()()2()()d i i nnx i x i i p x y x q x y x f x y x xS -=='=+-∑∑⎰记. (8.3.9)第2步 构造试探函数空间。
151第八章 常微分方程数值解在工程和科学技术的实际问题中,常常需求解常微分方程。
但由常微分方程理论可知,常微分方程中往往只有少数较简单和典型的方程可求出其解析解。
在大多数情况下,常微分方程只能用近似法求解。
这种近似解法可分为两大类:一类是近似解析法,如级数解法、逐次逼近法等;另一类则是数值解法,它给出方程在一些离散点上的近似解。
本章主要讨论一阶常微分方程的初值问题:()()⎪⎩⎪⎨⎧==0,y a y y x f dx dyb x a ≤≤ (8.1) 从理论上讲,只要方程中的()y x f ,连续且关于y 满足李普希兹(Lipschitz )条件,即存在常数L ,使()()2121,,y y L y x f y x f -≤-则常微分方程存在唯一解)(x y y =。
所谓微分方程数值解,就是求微分方程的解()x y 在一系列离散节点 b x x x x a n n =<<<<=-110处()i x y 的近似值i y ),,1,0(n i =. 相邻的两个节点之间的距离i i i x x h -=+1称为由i x 到1+i x 的步长,通常取为常数h 。
求数值解,首先应将微分方程离散化,常用的方法有: (1) 用差商代替微商 若用向前差商代替微商,即()()()()()i i i i i x y x f x y hx y x y ,1='≈-+ )1,,1,0(-=n i代入(8.1)中的微分方程,则得()1+i x y ()()()i i i x y x hf x y ,+≈152 记)(i x y 的近似值i y ,则由上式右端可计算出)(1+i x y 的近似值,即()i i i i y x hf y y ,1+=+ )1,,1,0(-=n i (8.2)(2) 数值积分法 利用数值积分法左矩形公式()()i i x y x y -+1=()()()i i x x y x hf dx x y x f i i,,1≈⎰+可得同样算法 ()i i i i y x hf y y ,1+=+(3) 用泰勒(Taylor )公式将函数)(x y 在i x 处展开,取一次Taylor 多项式近似,则得()()h x y x y i i +=+1()()i i x y h x y '+≈()()()i i i x y x hf x y ,+=从而也得到离散化得计算公式 ()i i i i y x hf y y ,1+=+§1 欧拉(Euler )方法1.1欧拉方法对一阶微分方程(8.1),把区间[]b a ,作n 等分:b x x x x a n n =<<<<=-110 , 则分点为 ih a x i +=, nab h -=),2,1(n i = 由以上讨论可知,无论用一阶向前差商,还是用数值积分法左矩形公式,或者用泰勒公式取前两项都可得到同样的离散化计算公式()i i i i y x hf y y ,1+=+并将初值条件代入,则得到数值算法:()()⎩⎨⎧=+=+a y y y x hf y y i i i i 01, ),2,1(n i = (8.3) 称其为欧拉方法。
几何上欧拉方法就是用一条折线近似表示曲线()x y y =(如图8-1)。
因此欧拉方法又称为欧拉折线方法。
153y图8-1 欧拉方法 1.2欧拉方法的误差估计定义1 假设)(i i x y y =为准确值,考虑计算一步所产生得误差,即用某种数值算法计算)(1+i x y 1+i y 所产生的误差()111+++-=i i i y x y R ,称为该数值算法的局部截断误差。
定义2 考虑用某种数值算法计算时,因前面的计算不准确而引起的准确解)(i x y 与数值解i y 的误差,()i i i y x y e -=称为该数值算法的整体截断误差。
设函数),(y x f 充分光滑,问题(8.1)的解()x y 在],[b a 上有二阶连续导数,由泰勒公式有()()h x y x y i i +=+1=()()()i i i y h x y h x y ξ''+'+221=()()i i i i y h y x hf y ξ''++221,所以154 ()111+++-=i i i y x y R =)(212i y h ξ'',1+<<i i i x x ξ (8.4) 定义3 如果一数值解法的局部截断误差为)(1+p h O ,则称该算法为p 阶算法。
当h 充分小时,由(8.4)知欧拉方法的局部截断误差为)(2h O ,因此欧拉方法是一个 一阶方法,计算结果的精度较差。
1.3 改进的欧拉方法由微分方程数值解的三种基本构造方法知,若取不同的差商(如向后差商),不同的数 值积分公式(如梯形公式),以及泰勒公式取前三项、四项等可得不同的算法。
如果用向后差商近似代替导数,则有()()()()()i i i i i x y x f x y hx y x y ,1='≈-- ),,1(n i =即 ()111,)()(++++≈i i i i y x hf x y x y )1,,1,0(-=n i所以有 ()111,++++=i i i i y x hf y y )1,,1,0(-=n i (8.5) (8.5)式称为隐式欧拉公式。
如果用梯形公式计算积分:()()()()()()[]11,,2,1+++≈⎰+i i i i x x x y x f x y x f hdx x y x f i i()()[]111,,2+++++=i i i i i i y x f y x f hy y (8.6) 且 ()111+++-=i i i y x y R = ()ξy h '''-3121(8.7)由于此方程为1+i y 的隐式方程,不易求解。
一般将其与欧拉方法联合使用,可得算法()()()()()()[]⎪⎩⎪⎨⎧++=+=+++++k i i i i i k i i i i i y x f y x f h y y y x hf y y 111101,,2, (8.8))1,2,1,0;,2,1,0(-==n i k按式(8.8)计算问题(8.1)的数值解时,如果每步只迭代一次,相当于将欧拉公式与梯形公式结合使用,即在实际计算中,当h 比较小时,常取一次迭代后的近似值()11+i y 为1+i y ,155于是有改进的欧拉方法⎪⎩⎪⎨⎧++=+=++++)]~,(),([2),(~1111i i i i i i i i i i y x f y x f h y y y x hf y y )1,2,1,0(-=n i (8.9) 例1 用欧拉方法和改进的欧拉方法求微分方程()[]7.0,0,10322∈⎪⎩⎪⎨⎧==x y xydx dy 的数值解(取h=0.1)。
解 由欧拉方法(8.3),得数值计算公式1+i y =i y +0.1×232i i y x 计算结果如表8-1由改进的欧拉方法(8.8),得数值计算公式⎪⎩⎪⎨⎧++=+=++++]~3232[05.0321.0~2112121i i i i i i i i i i y x y x y y y x y y 计算结果如表8-2表8-1x i 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y i 1.0000 1.0069 1.0208 1.0391 1.0628 1.0923 1.1269 1.1643 误差 0.0000 0.0037 0.0077 0.0993 0.0120 0.0151 0.0189 0.0222表8-2x i 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y i 1.0000 1.0033 1.0132 1.0292 1.0506 1.0773 1.1079 1.1422 误差 0.0000 0.0000 0.0000 0.0002 0.0020 0.0101 0.0122 0.0053例2 用欧拉法、改进欧拉法求微分方程数值解(h=0.1)。
156 ⎪⎩⎪⎨⎧=-='1)0(y y x y y解 由欧拉方法(8.3),得数值计算公式1+i y =i y +0.1⎪⎪⎭⎫⎝⎛-i i i y x y 由改进的欧拉方法(8.8),得数值计算公式⎪⎪⎩⎪⎪⎨⎧-+-+=-+=+++++)]~~()[(05.0)(1.0~11111i i i i i i i i i i i i i y xy y x y y y y x y y y计算结果如下§2 龙格-库塔(Runge -Kutta )方法2.1 泰勒展开法由于欧拉方法为一阶方法,为了提高算法的阶,有必要讨论更高阶的方法。
在泰勒展 开式中取更多的项,如取p +1项可得p 阶算法。
()p i p i i i i y p h y h y h y y !!221++''+'+=+157()()()i p p i yp h R ξ111!1++++= 其中()k i y )可用复合函数求导法则计算。
如p=2时得二阶泰勒方法()221(,),[]2!2i i i i i i i i i x y x y h h y y hy y y hf x y f f f +'''''=++=+++⋅2.2 龙格-库塔法为了避免计算高阶导数,龙格-库塔方法利用()y x f ,某些点处的值的线性组合构造计算公式,使其按泰勒公式展开后与初值问题解的泰勒展开式比较,有尽可能多的项相同。
龙格-库塔法的一般形式为:()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++==++++=+h y h x f K h y h x f K y x f K K a K a K a h y y m i m i m i i i m m i i μλμλ,,,222122111 (8.10) 下面以二阶龙格-库塔法为例说明龙格-库塔法的构造过程。
二阶龙格-库塔公式为(8.11)其中221,,λa a 为待定参数。
事实上,将K 2在()i i y x ,处按泰勒公式展开,则()i i i i i y x hf a y hK a hK a y y ,122111+=++=++()()()()()],,,,[2222h O y x f yy x hf y x f x hy x f h a i i i i i i i i +∂∂+∂∂+λλ ()()()⎪⎩⎪⎨⎧++==++=+1222122111,,hK y h x f K y x f K K a K a h y y i i i i i i λλ158 ()()()()()()322221,,,,h O y x f y y x f y x f x h a y x hf a a y i i i i i i i i i +⎥⎦⎤⎢⎣⎡∂∂+∂∂+++λ=()()()()()322221h O x y h a x y h a a x y i i i +''+'++λ 另一方面()()h x y x y i i +=+1=()()()()3221h O x y h x y h x y i i i +''+'+于是为使局部截断误差的阶尽可能高,应使⎪⎩⎪⎨⎧==+2112221λa a a 方程组有无穷多组解,取定参数则得到许多具体的二阶龙格-库塔公式。