1、航空航天中的计算方法航空航天中的计算方法授课教师：陈琪锋授课教师：陈琪锋中南大学航空航天学院中南大学航空航天学院第二部分第二部分 边值问题求解方法边值问题求解方法第第5章章 两点边值问题求解方法两点边值问题求解方法航空航天中的计算方法Page 42023-8-20内容提要内容提要5.1常微分方程边值问题的概念常微分方程边值问题的概念5.2打靶法打靶法5.3有限差分法有限差分法5.4有限元法有限元法1 Part 3:Two-Point Boundary Value Problems.2 David L.Darmofal,Computational Methods in Aerospace Eng
2、ineering(Lecture Notes),MIT,2005.Chap11,12.3 清华大学数学系编，现代应用数学手册清华大学数学系编，现代应用数学手册计算方法分册（计算方法分册（第十一章，常微分方程边值问题的数值方法），北京出版第十一章，常微分方程边值问题的数值方法），北京出版社，社，1990.航空航天中的计算方法Page 52023-8-205.1 常微分方程边值问题的概念常微分方程边值问题的概念 对于常微分方程：对于常微分方程：其中其中 ，x为标量，为标量，y和和 f 为为m维向量维向量。在。在上求解之需要上求解之需要m个定解条件，若定解条件的形式为：个定解条件，若定解条件的形式为
3、：其中其中 g为为m维向量。则该问题称为两点边值问题（维向量。则该问题称为两点边值问题（TPBVP，Two Point Boundary Value Problem）。）。如果边值条件形式可写为：如果边值条件形式可写为：其中其中gL和和gR的维数之和等于的维数之和等于m，则边界条件为分离的。，则边界条件为分离的。xxx()(,()yfy xydy d xa b,ab(),()0g yy 5.1 常微分方程边值问题的概念常微分方程边值问题的概念aa()0,()0LRgygy航空航天中的计算方法Page 62023-8-205.2 打靶法打靶法 以二阶系统为例，考虑边值问题：以二阶系统为例，考虑边
4、值问题：变换：变换：考虑初值问题：考虑初值问题：初值问题的解为：初值问题的解为：找到找到满足：满足：xxxxxa bab()(,(),(),(),()yfyyyAyB5.2 打靶法打靶法12yyyy xxxxxa bab1221211()()()(,),(),()yyyfyyyAyB xxa baa1221212,(,),(),()yyyfyyyAyxx12(;),(;)yyb1(;)yB 如何求如何求？航空航天中的计算方法Page 72023-8-20打靶法的几何解释：打靶法的几何解释：5.2 打靶法打靶法打靶：求解初值问题打靶：求解初值问题航空航天中的计算方法Page 82023-8-20
5、5.1.1 割线法割线法 以两个不同的以两个不同的值求解初值问题，得到两个解：值求解初值问题，得到两个解：根据初值条件知：根据初值条件知：假设假设 是是的线性函数，可取的线性函数，可取 为：为：迭代求解公式：迭代求解公式：结束条件：结束条件：5.2 打靶法打靶法 bbb1020101110(;)(;)(;)Byyy b1(;)yxx1011(;),(;)yyaa1011(;)(;)yyA bbb11111111(;)(;)(;)mmmmmmmByyy b111(;)myB 航空航天中的计算方法Page 92023-8-20割线法的几何解释：割线法的几何解释：5.2 打靶法打靶法线性近似：按割线
6、求根线性近似：按割线求根航空航天中的计算方法Page 102023-8-205.1.2 牛顿法牛顿法 求解非线性方程（组）：求解非线性方程（组）：在已知初值在已知初值0的处的处Taylor展开：展开：线性近似：线性近似：迭代求解公式：迭代求解公式：结束条件：结束条件：5.2 打靶法打靶法 bb110100(;)(;)yBy b1(;)yB bbb21111001010(;)(;)(;)yyyOB bb111(;)(;)mmmmyBy b111(;)myB b1(;)?y 航空航天中的计算方法Page 112023-8-20差分法求偏导数差分法求偏导数或采用其它数值微分方法。或采用其它数值微分方
7、法。f 可微时解偏导数微分方程可微时解偏导数微分方程微分方程对微分方程对求偏导：求偏导：5.2 打靶法打靶法bbb11101010(;)(;)(;)yyy xxa baa1221212,(,),(),()yyyfyyyAyxx12(;),(;)yy xa baa122121212,(;)0,(;)1yyyyyffyyyy初值问题，可解！初值问题，可解！（与割线法等价）（与割线法等价）割线代替切线割线代替切线航空航天中的计算方法Page 122023-8-20每一步迭代求解初值问题每一步迭代求解初值问题其中：其中：解得：解得：得到的终端值和对得到的终端值和对的偏导数：的偏导数：5.2 打靶法打靶
8、法axaaa1212122121212212,()(,),(),()0,()1yyyAyfyyyzzzffzzzzyy xa b,1212,yyzzxxxx1212(;),(;),(;),(;)yyzzbb11(;),(;)yy 航空航天中的计算方法Page 132023-8-20作业题作业题5：用牛顿打靶法求解两点边值问题用牛顿打靶法求解两点边值问题迭代初始条件取迭代初始条件取 。5.2 打靶法打靶法xxxxx2222sin(ln),12(1)1,(2)2yyyyy (1)0y 航空航天中的计算方法Page 142023-8-205.3 有限差分法有限差分法 以二阶系统为例，边值问题：以二阶
9、系统为例，边值问题：有限差分近似有限差分近似将区间将区间 等分为等分为N个子区间个子区间将将 在在xi处处Taylor展开：展开：xxxxxa bab()(,(),(),(),()yfyyyAyB5.3 有限差分法有限差分法 baa,0,1,2,ihxihiNN x()y xa b,2311()()()()()23!iiiiiiiy xy xy xxxyxxxyxxx用差分近似代替微分，将微用差分近似代替微分，将微分方程化为代数方程求解分方程化为代数方程求解航空航天中的计算方法Page 152023-8-20若取若取x=xi+1=x+ih：忽略二阶以上部分，得一阶导数的忽略二阶以上部分，得一阶
10、导数的前向差分前向差分近似：近似：若取若取x=xi-1=x-ih：忽略二阶以上部分，得一阶导数的忽略二阶以上部分，得一阶导数的后向差分后向差分近似：近似：5.3 有限差分法有限差分法1()()()iiiy xy xy xh 23111()()()()()26iiiiiy xy xy x hh yxh yx 23111()()()()()26iiiiiy xy xy x hh yxh yx 1()()()iiiy xy xy xh 一阶精度一阶精度一阶精度一阶精度航空航天中的计算方法Page 162023-8-20 xi+1和和xi-1在在xi处的处的Taylor展开相减，忽略三阶以上部分，得展
11、开相减，忽略三阶以上部分，得一阶导数一阶导数的中心差分的中心差分近似：近似：xi+1和和xi-1在在xi处的处的Taylor展开相加，忽略四阶以上部分，得展开相加，忽略四阶以上部分，得二阶导数二阶导数中心差分中心差分近似：近似：三阶导数的中心差分近似？三阶导数的中心差分近似？5.3 有限差分法有限差分法11()()()2iiiy xy xy xh 二阶精度二阶精度112()()2()()iiiiy xy xy xyxh 二阶精度二阶精度航空航天中的计算方法Page 172023-8-20 xi+1和和xi-1在在xi处的处的Taylor展开相减，忽略五阶以上部分：展开相减，忽略五阶以上部分：x
12、i+2和和xi-2在在xi处的处的Taylor展开相减，忽略五阶以上部分：展开相减，忽略五阶以上部分：三阶导数的中心差分近似：三阶导数的中心差分近似：四阶导数的中心差分近似：四阶导数的中心差分近似：5.3 有限差分法有限差分法21123()2()2()()()iiiiiy xy xy xy xyxh 二阶精度二阶精度35111()()2()()()3iiiiy xy xy x hh yxO h35228()()4()()()3iiiiy xy xhy xh yxO h(4)21124()4()6()4()()()iiiiiiy xy xy xy xy xyxh 二阶精度二阶精度航空航天中的计算
13、方法Page 182023-8-20有限差分法解微分方程两点边值问题有限差分法解微分方程两点边值问题微分方程微分方程离散化，将区间离散化，将区间 等分为等分为N个子区间：个子区间：在节点上应用中心差分公式，得到代数方程组：在节点上应用中心差分公式，得到代数方程组：xxxxxa bab()(,(),(),(),()yfyyyAyB5.3 有限差分法有限差分法 x1111202(,),1,2,12,iiiiiiiNyyyyyfyiNhhyAyB baa,0,1,2,ihxihiNN xa b,航空航天中的计算方法Page 192023-8-20有限差分法解微分方程两点有限差分法解微分方程两点边值问
14、题的几何解释边值问题的几何解释5.3 有限差分法有限差分法离散点：微分用有限差分近似离散点：微分用有限差分近似航空航天中的计算方法Page 202023-8-20例例5.1：用有限差分法求解两点边值问题：用有限差分法求解两点边值问题取离散化区间取离散化区间h=0.1，N=10。xxxxx2222sin(ln),12(1)1,(2)2yyyyy 5.3 有限差分法有限差分法 x1111202(,),1,2,12,iiiiiiiNyyyyyfyiNhhyAyB xxxxx111122202sin(ln)2221,21,1,2,1iiiiiiiiiiNiyyyyyyhhyyihiN 航空航天中的计算
15、方法Page 212023-8-20线性方程组：线性方程组：即：即：5.3 有限差分法有限差分法 xxxxxx2211220sin(ln)12 111,21,1,2,1iiiiiiiiNihhhhyyyyyihiNxxxxxxxx10112222111()1()sin(ln),()()1NNNNHYGhgyyyghYGgyhgy 航空航天中的计算方法Page 222023-8-205.4 有限元法有限元法 以二阶系统为例，考虑边值问题：以二阶系统为例，考虑边值问题：5.4.1 投影类方法的基本思想投影类方法的基本思想 以一简单函数以一简单函数 近似近似y(x)，给出连续近似解，例如：，给出连续
16、近似解，例如：一般形式：一般形式：，已知，已知，待定。待定。残差：残差：某种意义上使残差最小，则得到某种准则下最佳的近似解。某种意义上使残差最小，则得到某种准则下最佳的近似解。xxxxxa bab()(,(),(),(),()yfyyyAyB5.4 有限元法有限元法x()yxxxxx()()(,(),()Ryfyy x10()Nkkkyc x x0()kkkyc x x 10()()Nkkkyx k()kx航空航天中的计算方法Page 232023-8-20区间残差平方和最小：最小二乘法区间残差平方和最小：最小二乘法若干特定点处残差为零：配点法若干特定点处残差为零：配点法加权残差为零：加权残差
17、法加权残差为零：加权残差法Galerkin法：法：。5.4 有限元法有限元法x012,min()NbaRdx 12()0,1,2,iNRiNab xx()()0,1,2,biawRdxiN ()iw x计算复杂，不常用计算复杂，不常用为权函数为权函数()()iiw xx 航空航天中的计算方法Page 242023-8-20例例5.2：考虑两点边值问题：考虑两点边值问题解析解为：解析解为：试用配点法和加强残差法求解该问题近似解。试用配点法和加强残差法求解该问题近似解。x50,11(1)100,(1)100 xyeyy 5050sinh(1),10050cosh(1)11sinh(),cosh()
18、22xxxxxyeaxbabxeexee 5.4 有限元法有限元法航空航天中的计算方法Page 252023-8-205.4 有限元法有限元法解析解解析解航空航天中的计算方法Page 262023-8-20设近似解的形式：设近似解的形式：基函数的选择示例：基函数的选择示例：为满足边值条件要求为满足边值条件要求取二次函数取二次函数以及三次项以及三次项N=2（1）配点法）配点法近似解的残差近似解的残差令令N个点处残差为零求解系数，如个点处残差为零求解系数，如x 1()100()Niiiyx 5.4 有限元法有限元法(1)0i1()(1)(1)xxx2()(1)(1)xxxx线性函数不满足线性函数不
19、满足xx12()502650 xxRyee x1 3 121125cosh,25sinh33配点？配点？航空航天中的计算方法Page 272023-8-20（2）加权残差法）加权残差法要求：要求：Galerkin法，取法，取即：即：5.4 有限元法有限元法111212587575,22eee1122()(),()()wxxwxx xx()()0,1,2biawRdxi x1121(1)(1)(2650)0 xxxedx x1121(1)(1)(2650)0 xxxxedx 航空航天中的计算方法Page 282023-8-205.4 有限元法有限元法-1-0.8-0.6-0.4-0.200.20
20、.40.60.8195100105110115120125130 exactCollocationGalerkin配点法、配点法、Galerkin加权残差法与精确解的比较加权残差法与精确解的比较航空航天中的计算方法Page 292023-8-205.4.2 有限元法的基本思想有限元法的基本思想 将区域（区间）划分为小的单元，在单元上表示近似解将区域（区间）划分为小的单元，在单元上表示近似解以及求残差加权积分。以及求残差加权积分。第第i个单元，个单元，Ei，在每个单元上解用多项式近似；在每个单元上解用多项式近似；在每个单元上计算加权残差；在每个单元上计算加权残差；根据各单元满足的方程确定多项式近
21、似解的系数。根据各单元满足的方程确定多项式近似解的系数。5.4 有限元法有限元法xxx1ii xxx1iii 局部近似，分段光滑局部近似，分段光滑可以用简单的低阶近似可以用简单的低阶近似航空航天中的计算方法Page 302023-8-205.4.3 有限元法：线性元为例有限元法：线性元为例 解在每个单元上采用解在每个单元上采用x的线性函数近似表示。的线性函数近似表示。5.4 有限元法有限元法航空航天中的计算方法Page 312023-8-20 线性函数具有线性函数具有2个自由度：由两个端点的函数值确定；个自由度：由两个端点的函数值确定；N个线性单元，近似连续函数，个线性单元，近似连续函数，N+
22、1个自由度：由个自由度：由N+1个节点的函数值唯一确定。个节点的函数值唯一确定。设近似解表达为：设近似解表达为：由由 可知可知5.4 有限元法有限元法x 11()()Niiiyx xx 11()(),1,2,1NjiijiyjN x121,1,2,1(),1,2,1iNjjL yyyiNyyjN 线性函数线性函数航空航天中的计算方法Page 322023-8-20 近似解可以用节点基函数表示为：近似解可以用节点基函数表示为：节点基函数节点基函数 在节点在节点 j 处：处：由于由于 近似表达式中近似表达式中 取值的任意性，可知：取值的任意性，可知：5.4 有限元法有限元法xx11()()Niii
23、yy xxxx111111()()()()NjiijijNiijjjjiijiijyyyyy 节点基函数的特性节点基函数的特性x()yx1,()0,ijijij 航空航天中的计算方法Page 332023-8-20 对于线性元，节点基函数在每个单元内是线性函数：对于线性元，节点基函数在每个单元内是线性函数：5.4 有限元法有限元法xxxxxxxxxxxxxxxxx11111110()()()0(2,3,)iiiiiiiiiiiiN 第第i 个节点基函个节点基函数的几何表示数的几何表示xxxxxx111iiiiii xxxxxxxxx211212()()0 xxxxxxxxx110()NNNNN
24、N 航空航天中的计算方法Page 342023-8-20节点函数值的求解：加权残差节点函数值的求解：加权残差Galerkin法法 近似解的节点基函数表示：近似解的节点基函数表示：Galerkin法求法求 ，即解方程组：，即解方程组：其中：其中：（对于示例二阶系统）（对于示例二阶系统）即：即：分部积分：分部积分：5.4 有限元法有限元法xxx11()()()Niiiyy x()iyxx()()0,1,2,1biaRdxiN xxxxx()()(,(),()Ryfyy ()0,1,2,1bbbiiiaaayydxfdxiN 0,1,2,1bbiiaadyfdxiN 航空航天中的计算方法Page 3
25、52023-8-20在第在第i个个单元单元内内，Ei：5.4 有限元法有限元法单元内单元内 为连续函数为连续函数x()yxxxxxxxxxxxxx1111101()10(2,3,)iiiiiiiiiiN xxxxxxx112121()0 xxxxxxx110()1NNNNN xxxxxx11111()()()()()Niiiiiiiiiiyyyyyy xxx1ii 航空航天中的计算方法Page 362023-8-20方程组：方程组：当当 时方程为：时方程为：当当 时时方程为：方程为：当当 时时方程为：方程为：5.4 有限元法有限元法1i xxx211122()0 xaxyyyfdx ()0,1
26、,2,1bbbiiiaaayydxfdxiN 2,3,iN 11110iiiixxiixxydxfdx 2211110 xxaxxyydxfdx xxxxxxxxxx111111111()()0iiiiiiiiiiixxiiiixxyyyfdxfdx 1iN 11110NNNNxxbNNxxyydxfdx xxx11()0NNxNbNNNxyyyfdx 航空航天中的计算方法Page 372023-8-20 给定边值条件：给定边值条件：上述上述N+1个方程可解出个方程可解出N+1个未知量：个未知量：最终得到问题的解：最终得到问题的解：5.4 有限元法有限元法xx11(),()NyAyB xaxb
27、11,N xxx11()()()Niiiyy xxx23,(),(),(),aNbyyyyy航空航天中的计算方法Page 382023-8-20例例5.2（续）：两点边值问题（续）：两点边值问题用线性元和加权残差用线性元和加权残差Galerkin法得到法得到N+1个方程：个方程：给定：给定：，可求解剩余，可求解剩余N+1个未知数。个未知数。x50,11(1)100,(1)100 xyeyy 5.4 有限元法有限元法 xxxxxxxxxxxxxxxx2111111221111111150()50()50(),2,3,50()0iiiiNNxxaxiiiiiiixxxxiiiixxxxNbNNNx
28、yyye dxyyye dxe dxiNyyye dx xxx(,(),()50 xfyye 11100,100Nyy 航空航天中的计算方法Page 392023-8-20高斯求积高斯求积 将积分表示为被积函数在若干点处的函数值加权和：将积分表示为被积函数在若干点处的函数值加权和：若适当选取若适当选取 和和 ，可使公式对次数，可使公式对次数 2N+1的多项式被积的多项式被积函数均精确成立（具有函数均精确成立（具有2N+1次代数精度），这类公式称为次代数精度），这类公式称为高斯求积公式。节点高斯求积公式。节点 称为高斯点。称为高斯点。5.4 有限元法有限元法110()()Niiigdg iii航
29、空航天中的计算方法Page 402023-8-20例：例：N=0时，只能实现对时，只能实现对1次多项式精确满足，对比系数，得：次多项式精确满足，对比系数，得：N=1时，不大于时，不大于3次的多项式可精确次的多项式可精确满足，对比系数，得：满足，对比系数，得：N=，高斯点和系数可查数学手册，高斯点和系数可查数学手册。5.4 有限元法有限元法0020 1001()()gdg 010111313 100111()()()gdgg 110()()Niiigdg 注意积分区间变换到注意积分区间变换到-1,1！航空航天中的计算方法Page 412023-8-205.4.4 有限元法拓展有限元法拓展高次元：
30、分段高次插值高次元：分段高次插值多维问题多维问题5.4 有限元法有限元法航空航天中的计算方法Page 422023-8-20作业题作业题6：用有限元法，采用线性元，求解两点边值问题：用有限元法，采用线性元，求解两点边值问题：将区间等分为将区间等分为10个单元，编写个单元，编写MATLAB程序求解。程序求解。5.4 有限元法有限元法x50,11(1)100,(1)100 xyeyy x40sin0,00.5(0)(0.5)0yyyy 航空航天中的计算方法Page 43人有了知识，就会具备各种分析能力，人有了知识，就会具备各种分析能力，明辨是非的能力。明辨是非的能力。所以我们要勤恳读书，广泛阅读，所以我们要勤恳读书，广泛阅读，古人说古人说“书中自有黄金屋。书中自有黄金屋。”通过阅读科技书籍，我们能丰富知识，通过阅读科技书籍，我们能丰富知识，培养逻辑思维能力；培养逻辑思维能力；通过阅读文学作品，我们能提高文学鉴赏水平，通过阅读文学作品，我们能提高文学鉴赏水平，培养文学情趣；培养文学情趣；通过阅读报刊，我们能增长见识，扩大自己的知识面。通过阅读报刊，我们能增长见识，扩大自己的知识面。有许多书籍还能培养我们的道德情操，有许多书籍还能培养我们的道德情操，给我们巨大的精神力量，给我们巨大的精神力量，鼓舞我们前进鼓舞我们前进。
