第二章发展方程的有限元分析

第二章发展方程的有限元分析
W. B. J. ZIMMERMAN，B. N. HEWAKANDAMBY
Department of Chemical and Process Engineering, University of Sheffield,
Newcastle Street, Sheffield S1 3JD United Kingdom
E-mail:
科学研究和工程应用中的偏微分方程（PDE）多源自复杂的平衡方程。

常见的偏微分方程主要来自质量守恒、动量守恒、组分守恒和能量守恒定律。

由于这些守恒定律是整个域上的积分方程，所以在连续性假设下，偏微分方程很容易用有限元方法近似描述。

本章介绍了COMSOL Multiphysics中典型的三种不同类型“时间－空间”系统偏微分方程——椭圆方程，抛物线方程和双曲线方程。

本章还对有限元方法进行了总体介绍，结合应用实例讲解有限元方法精确计算的特性，更深层次的内容将在后续章节中引出。

1. 简介
在科学研究和工程应用中常会遇到满足守恒定律约束的偏微分方程，通常以积分形式出现。

所有由质量守恒、动量守恒、组分守恒和能量守恒控制的传递现象都会产生连续逼近的偏微分方程。

相信化工人员对传热、传质和动量传递现象不会感到陌生。

与前一章COMSOL Multiphysics化工实例中介绍的零维、一维空间系统相比，化工课程中通常不会出现超过二维或三维的偏微分方程计算。

从文献[1]中找到一个非常珍贵的例子。

实际上，很多常见的化工模型和公式都是由实际过程中更高维数的动力学过程简化而来的。

流体动力学中的阻力系数，传热传质系数，多相催化的Thiele模型，精馏塔设计中的McCabe-Thiele图等许多描述高维数系统传递现象或非稳态动力学过程的技术，都是半经验性的方法，也许可以用偏微分方程来描述这些过程，但是由于基本物理、化学现象的复杂性，这些方程通常很难求解。

所以对于初步的设计计算，这些快速计算的简化方法很受欢迎，但是对于细节设计、设计翻新、过程分析和优化过程，只有简化方法是不够的。

在基础科学研究过程中，这些方法仍然不断地从化学工程师传给生物学家、材料学家等，逐渐应用到各个领域的工作中。

但是计算流体动力学（CFD）的出现彻底改变了这些方法在传导模型中的地位。

虽然CFD在传导现象可视化和量化方面具有独特优势，但是唯象方法对于描述分布系统模型仍然有非常重要的作用。

COMSOL Multiphysics不是一个“商业CFD软件”，但是也可以做一些CFD 计算。

它包含一些通用的CFD软件包，在某些模型支持上具有其独特的优点。

对于CFD方法，大多数过程工程师希望有对湍流和燃烧模型的支持。

但是COMSOL Multiphysics不同，擅长多物理场计算。

除了CFD的传统传递现象，COMSOL Multiphysics还包含了电动力学、磁动力学、结构力学等实例模式，可以对这些现象进行模拟计算。

实际上COMSOL Multiphysics的最大优点在于——首先，用户自定义编程容易，可以轻易建立用户自定义的模型，通过改变变量系数、边界条件、初始条件，还可以同时耦合多个物理场（甚至不在同一个域上的物理场）；其次，基于MATLAB（或COMSOL Script），可以实现对复杂模型模拟的所有编程功能，可以把COMSOL Multiphysics看作是一个方便的高阶有限元编程和分析软件。

上一章中，我们介绍了一些用户自定义编程分析的功能，本章将介绍COMSOL Multiphysics在更高维数偏微分方程系统有限元建模上的核心
优势。

多物理场、扩展多物理场和无PDE 约束的处理将留给后续章节。

2． 偏微分方程
偏微分方程可以根据阶数、边界条件类型和线性度（线性、非线性和准线性）进行分类。

令人惊讶的是，大多数科学研究和工程应用中碰到的偏微分方程都是二阶的，即最高导数项是二阶偏导数。

这是偶然的么？这一巧合往往是在各自领域分别解释的。

但是，最近Frieden[2]证明所有已知物理定律都可以由最小Fisher 信息原理推导出，而该原律通常将一个二阶项作为最高项引入相应的物理定律中——包括从Schrodinger 波函数方程到经典电动力学等各种物理定律。

所以二维和三维的二阶“空间－时间”系统求解和分类在科学研究和工程中有着广泛的应用，这非常重要！出于这个原因，且有限元方法本质上非常适合处理二阶系统，所以有限元是有着广泛应用的技术。

本章我们将集中介绍二维和三维空间的二阶系统。

这里有三个典型的例子，在课本中经常见到，它们是：
222222
2222Laplace (): 0 (): (): u u x y
u u t x u u t x
∂∂+=∂∂∂∂=∂∂∂∂=∂∂方程椭圆 扩散方程抛物线 波方程双曲线 下面用COMSOL Multiphysics 对这几种偏微分方程进行建模，同时在每个例子中我们将引入一个新的或特殊的COMSOL Multiphysics 功能。

这些椭圆、抛物线和双曲线方程可以归结为一种通用的线性二阶偏微分方程，具有两个独立变量和一个非独立变量：
2222220u u u u u a b c d e fu g x x y y x y ∂∂∂∂∂++++++=∂∂∂∂∂∂ (1)
方程中的系数只是独立变量x 和y 的函数，或者是常数。

根据以下原则分为三种经典形式：
222 : 0
: 0: 0
b a
c b ac b ac -<-=->椭圆抛物线双曲线
这种分类方式大致反应了域中的信息流。

例如，对于椭圆方程，边界点的信息会快速地传到所有内部点。

所以椭圆方程是“非本地的”，意味着很远处的信息会对给定位置有影响，而“本地的”则表示只有附近信息才对其有影响。

抛物线系统“扩散”信息，即信息会向所有方向传播。

双曲系统“传播”信息，也就是说已经接收到信息的区域和将要接受到信息的区域有一个明显的划分。

如果系统是线性的或者准线性的（即一些系数依赖于非独立变量，或者乘以一个更低阶的偏导数），这种分类方式和信息传播方式可以作为对二阶系统的一个引导。

但是在非线性系统中，非线性破坏了信息传递结构，信息可能发生“异常”，也就是说，不再传递、超出给定范围、来自噪声或者在给定窗口内因丢失初始条件而消失。

2.1 泊松方程：椭圆偏微分方程
拉普拉斯方程经过适当变化就成为泊松方程：
2()
∇=(2)
u f x
我们在第一章描述非均匀介质中具有分布热源的热传导问题时见过该方程的一维形式，但是这里的热传导率是均匀的。

下面给出(2)式的另一个侧面，应该注意到，该方程满足给定涡量曲线的流函数方程：
2()x
ψω
∇=-(3) 有两种较为容易的常见涡量类型——兰金涡（某一区域内的涡量是常数）和点源涡（涡量快速降低，可以理想为一个点涡）。

你可能会对这两种类型涡产生的流线感兴趣。

下面我们来研究一下这些流线。

打开COMSOL Multiphysics模型导航栏。

按照表1的步骤求解具有单位分布源项（常涡量）的泊松方程。

该模型具有一个独立变量u，二维空间坐标x和y。

我们绘制一个实心圆作为研究域，COMSOL Multiphysics将边界分为四部分，默认网格划分762个。

等流函数线如图1所示。

图1 光滑圆形域中的常涡量流线
流线用等高线的形式显示，图1通过“Menu”菜单中“Export：Image”选项导出。

显然流线是一组同心圆。

整个边界也是条流线（ψ＝0），最大流函数出现在圆心（ψ＝0.250），所以由常涡量引入的体积流量为0.25。

不改变几何形状，将
网格数加密到3045个单元。

由于流函数等高线相减的差值就是流线间的体积流率，流线的切线方向就是流动方向，所以图1的含义就很清楚了——常涡量导致圆形域中的环形流动。

涡量从边界扩散到每个地方，导致边界上产生最大流速，圆心速度最小。

下面来看域中的点源情况。

我们将在域中画一个非常小的圆形作为第二个域，来近似点源情况。

在第二个域中设定f=1/πR2，这里R是小圆的半经，小圆外侧f＝0。

这样f积分后为单位1，当R→0时，f趋近于Dirac delta函数。

但是极限将通过细化网格和使用弱形式来逼近。

打开“Draw”下拉菜单，在域中定义一个点。

在域中明确绘制一个容易区分的点，会出现两种情况——可能（但不是必须）会将点贡献或点约束施加在方程系统上，也可能点本身成为产生网格的节点，约束网格。

在网格节点容易（可能计算的更准确）求解独立变量的值，因为不需要进行值的网格内插。

在该点附近可能需要更精细（或粗糙）的网格划分。

按照表2的步骤来实现。

Draw菜单Specify objects：Point
设置点x=0, y=0，完成
Physics：Point settings Point selection：设置为3（原点）在week term项输入u_test
完成
Physics菜单：Subdomain settings 选择域1
设置c=1；f＝1 完成
Meshing 点击三角形按钮绘制网格
Solve菜单：Solver parameters Advanced选项卡：改变解的形式为weak 完成
Post-processing 点击原点。

报告栏显示：Value：0.805017，表达式：u，
位置：(0,0)
图2位于圆形域圆点处的点涡量产生的流线
图3位于圆形域圆点处的点涡量产生的流线
图2给出了表2中定义的弱的点源项产生的流线。

从域中流函数的值能够发现，体积流率要比均匀分布源项情况下大很多。

从本质上讲这是一组围绕圆点的等高线，不太满足圆心对称，当然在趋近圆心处同心圆消失了。

读者可能非常想知道“u_test ”是什么。

如何能够将“弱的点条件”转化为点源。

关于弱条件和弱约束的讨论我们将留在本章后面，等我们更加详细的介绍完有限元方法以后再加以介绍。

但是这里提一下，它需要满足有限元方法的一个强项，就是不连续函数，如Dirac delta 函数、步进函数的半分析处理。

这些问题起源于工程动力学，例如合成材料或普通材料的性能在表面会发生变化。

将网格细化到2986个单元，并不能使圆形等高线看起来更为圆滑，但是最大流函数值从0.805增大到0.915。

这是由细化圆心位置网格导致的。

打开“Mesh ”菜单，选择“Parameters ”选项。

在“Point ”选项卡中输入顶点3的最大基元尺寸为0.001。

重新划分网格，会产生1658个基元，但是原点处的网格更密，计算得到圆点处ψ＝1.57。

图2中最大流函数为1.57。

将网格增加到6632，最大流函数为1.68。

类似增大网格求解，26528个网格对应流函数1.79，106112个网格对应流函数1.90。

虽然不知道网格是否能够收敛，但是由于涡量随着离源点距离的增加快速降低，本质上流线分布已经收敛。

练习2.1
求解当涡量随半经指数降低时的流线，即： 220exp(x y ωω=-+
2.2 扩散方程：抛物线偏微分方程
一维非稳态扩散方程的形式能够很好地适用于以下三种常见的传递性现象：
222222 p c c D t x
T k T t c x
t x
ρωων∂∂=∂∂∂∂=∂∂∂∂=∂∂ 质量扩散： 热传导：二维涡量输运： 这里c ，T 和ω分别是浓度、温度和二维流动中z 方向的涡量，它们相应的扩散系数为D ，α和ν。

该方程已经在本科课程中详细介绍过。

它可以通过傅立叶和拉普拉斯变换求解，初始和边界条件下的相似解是否失效取决于变量Dt x =η。

这里没有有限元方法求解的机会——只有求解这类老问题的老办法，对吗？不对。

COMSOL Multiphysics 仍然可以改进、计算这类问题，这一点常被忽略。

COMSOL Multiphysics 非常适用于非常数系数情况，即传递性能取决于场变量。

例如，对于低压和高温情况，气体必须满足理想气体定律：
nM PM V RT
ρ=＝ (4) 这里R 是气体常数，M 是组分相对原子质量。

在这种情况下，很少有气体热容能够保持常数。

例如，在一定温度范围内，CO 2气体的热容可以用温度的二次函数近似(MJ/kg mol ℃)：
5236.110.04233 2.88710p c T T -=+-⨯
(5)
需要满足： 372()36.11(273)()1 1.172107.99510p k kR f T c PM
T f T T T ρ--=+=
+⨯-⨯ (6)
对于时间和位置满足： 2
36.11kR
x t X PML L τ== (7)
带入热传导方程得到简化形式为： 22()T T f T X τ∂∂=∂∂ (8)
下面来考虑在CO 2气体层中的热传导，设气体层长度为L ，水平边界保持恒温400℃和500℃。

打开COMSOL Multiphysics模型导航栏，按照表3的步骤建立模型。

该模型给出一个独立变量u和一维空间坐标x。

点击工具栏上网格划分按钮，然后再次加密网格两次，使网格单元达到60个。

点击“＝”按钮进行求解（稳态线性求解器）。

稳态解不随初始条件发生变化，因为经过足够长的时间后，累计项逐渐消失。

当然在稳态求解器下是这样。

所以我们可以用依赖于温度的扩散系数来研究对瞬态解的影响。

所以改变初始条件为u(t0)＝400℃，即t＝0时左侧边界条件跃升至u＝500℃。

打开“Solver”菜单选择“Parameters”选项，会弹出求解器变量设置对话框。

按表4内容进行瞬态分析设置。

经过后处理生成图4和图5。

图4给出了不断逼近稳态的温度曲线（几乎线性），图5给出了整个区域中间点的温度变化。

求解（＝）
Post-processing：
接受默认设定，完成
Domain Plot
Point选项卡，设定x：0.5，完成Post-processing：
Cross-section plot
parameters Array
图4 从t＝0到t＝0.2的温度曲线
图5 x ＝0.5处的温度变化
图4和图5给出了向前扩散的速率，问题是非线性介质在其中起什么作用。

特别是由于扩散率随温度升高而增大，我们发现与常数扩散率相比，曲线能够更快达到稳态。

由于表4中没有出现自相似项Dt
x =η，所以高温比低温更快达到稳态线性曲线。

图5给出了域中心点处温度升至稳态值的过程，曲线呈S 形，比预计要快。

练习2.2
设常数扩散率f (T )=1求解图4和图5，这些曲线有什么区别？
2.3 波动方程：双曲偏微分方程
一维正则波动方程已经研究的非常透彻，尤其在化工实例中经常见到。

最典型的就是化工课程中容易忽视的声波。

这是否意味着化学和过程工业中波动方程不重要呢？实际上总是会遇到很多复杂的波动方程或者非线性发展方程。

例如，前沿研究经常关注的化学波反应器，冷凝器表面的波，旋转喷油嘴，分裂蒸馏塔质量传递影响，声学，超声波和声化学等等。

但是碰巧化学工程师很少学习波的知识，市面上也很难找到经典的课本，分析声学在化工单元中的应用。

我们将一维正则波动方程留作练习，这里以Korteweg-de Vries 波传播方程为例加以介绍：
3360u u u u t x x ∂∂∂++=∂∂∂ (9) 这是薄层表面波的波方程——常见的用石头打水飘就是这样一个过程——在物
理学中广泛存在，因为它可以在任何波传播介质中产生，不会消散，从某种意义上讲是非常“微弱的”。

作者最近在地球物理[3]和表面涂层应用[4]方面对这类波方程加以推导。

在这两方面工作中，都使用特殊的线性化方法(MOL)模拟作为数值分析的指导。

线性化方法使用有限差分方法离散偏微分项，所以在网格点产生了一套可以同时求解u 值的常微分方程。

COMSOL Multiphysics 用有限元方法求解时间依赖的偏微分方程系统，本章后续章节还会遇到类似情况。

与使用常微分方程计算网格点的u 值不同，从COMSOL Multiphysics 求解常微分方程中产生的“自由度”，更容易得到任何位置独立变量u 的值。

出于技术原因，我们不考虑用FEMLAB 建模计算从[3][4]中推导出的类似与(9)的波方程——我们希望严格控制数值积分过程，使用刚性系统或者全隐式进行时间积分，因为这类方程很容易导致数值非稳定性。

但是COMSOL Multiphysics 具有强有效的时间积分求解器，会自动选择刚性求解器。

所以仅仅出于好奇，我们决定用COMSOL Multiphysics 求解孤波传递的Korteweg-de Vries 方程。

我们选择COMSOL Multiphysics 而不用FEMLAB 的一个重要原因在于对高
阶偏导数项的处理。

据我所知，FEMLAB 不能处理类似33x
u ∂∂的项。

在COMSOL Multiphysics 中也没有定义u xxx ，但是定义了u xxt ，使得我们可以通过一些办法来构造u xxx 。

我的办法是引入一个辅助变量v ，生成一个与(9)等价的系统：
222[3]0u u v t x u v x ∂∂++=∂∂∂=∂ (10)
第一个方程是通用守恒形式，第二个方程是一个关于v 的椭圆方程。

所以整个系统是一个“微分－几何”方程，不是所有的变量都需要“发展”。

但是初始系统具有双曲特性。

由于假设波传播没有方向性，理想情况下应当在无限介质中模拟。

因为没有哪台计算机有无限的内存，所以提供足够大介质使得边界效应可以忽略的情况就显得更有意义，于是引入周期性边界条件。

这里我们将根据方程(9)模拟一个孤波的传播。

如果这些孤立或者本地的结构足够紧凑，域足够长，波就不会和它的周期性反射波相互作用。

通过结果分析我们将看到这一点。

方程(9)的孤波解析解为：
223(,)2sech (4)u x t a ax a t =-
这里a 是振幅。

注意到传播速率（系数x 与t 的比值）依赖于a 2，因此波振幅越大传播速度越快。

这是由Korteweg-de Vries 方程的非线性引起的。

有必要知道初始条件并进行模拟，用二阶偏微分项定义v ：
44(,0)4(cosh(2)2)sech ()v x t a ax ax ==-
打开COMSOL Multiphysics 模型导航栏。

表5中包含建立孤波传递模型的具体步骤。

表5 非线性波方程举例——Korteweg-de Vries 孤波传递
图6 KdV方程孤波传递解的初始波形
图7 KdV方程孤波传递解拉普拉斯算子v的初始波形
图6和图7分别给出了孤波解和方程(10)的拉普拉斯算子的初始波形。

表5
建立了一个周期性边界条件。

应当注意到，对于Korteweg-de Vries 系统(9)，正振幅的波向右移动，负振幅的波向左移动。

由于这里只有时间的一阶偏导数项，所以Korteweg-de Vries 方程不可能使波向正则波动方程那样同时向两个方向传递。

练习2.3 正则波动方程，2222x
u t u ∂∂=∂∂，具有相应的经典偏微分方程应用模式。

该应用模式具有二阶空间和时间偏导数。

改变边界条件（Neumann 边界条件），对于正则波动方程建立表5中的周期性边界条件。

尝试初始条件u (t 0)=sin(20*Pi *x )和u _t (t 0)=-10*Pi *cos(10*Pi *x )。

你预计u (t )和u _t (t )的结果如何？计算结果是否与你预料的不同？注意COMSOL Multiphysics 有预定义的pi 。

练习2.4
当a ＝2时，求解Korteweg-de Vries 模型。

理论上讲，振幅因子会导致传播速度增大a 2倍。

模拟结果如何？对克服这个问题有什么建议？
3．有限元方法
到目前为止，你一定很想知道COMSOL Multiphysics 是如何求解偏微分方程系统的。

有限元分析已经存在了几十年，从20世纪80年代开始就出现了相应的商业软件，Reddy[5]的书对此给出了很好的介绍。

这里不是想详细介绍有限元方法，也不是仔细讲解COMSOL Multiphysics 求解是如何实现的，而是希望对有限元方法计算类型有一个大致印象，理解为什么说COMSOL Multiphysics 是实现有限元方法的一个很方便的工具。

有限元方法的本质是将任何场变量约束以弱形式来表示。

理解什么是弱形式（以及为什么数学上要求必须是弱形式），就必须要理解系统约束的强形式是偏微分方程系统和适当的边界条件。

它为什么强？由于场变量需要连续并且要有与方程阶数相同的连续偏导数。

这就是强约束。

而弱形式对函数满足的约束条件限制要相对弱一些——不连续但是可积。

来看一个偏微分方程和它的等价弱形式。

考虑在三维空间域Ω中的具有一个独立变量u 的稳态偏微分方程，通用形式为：
()()u F u ∇⋅Γ= (11)
我们假设测试函数v 是域Ω中定义的一个任意函数，并且满足函数v ∈V 。

方程
(11)两侧同乘以v ，并在整个域内积分有：
()()v u dx vF u dx ΩΩ∇⋅Γ=⎰
⎰ (12) 这里dx 表示体积微元。

根据散度定律，我们得到： ()()v ndS v u dx vF u dx ∂ΩΩΩΓ⋅-∇⋅Γ=⎰⎰⎰ (13)
当偏微分方程被Neumann 边界条件约束时，方程(13)中的边界项消失。

这也是有限元方法有Neumann 自然边界条件的一个原因。

回顾第一章，我们发现有限微分法有自然的Dirichlet 边界条件。

这导致体积积分变为：
[()()]0vF u v u dx Ω-∇⋅Γ=⎰ (14)
这必须对所有的v ∈V 都满足。

现在来看有限元和基底函数。

假设u 可以分解为一系列基底函数：
()()i i i u x u x φ=∑
(15)
例如，如果i φ是正弦、余弦等基本渐进谐波函数，那么(15)就是由一系列傅立叶函数构成的。

相反，在有限元方法中，要选择那些只在某一基元中作用的基底函数，即除了某一基元，在其它基元中都为零。

图8一维空间临近基元的两段线性基底函数
图8给出了一维空间中两个拉格朗日线性基底函数的例子。

显然，任意函数u (x )都可以由分段线性基底函数和足够小的基元进行任意精确度的近似。

基底函数可以使用更高阶数，但是相应的每个基元也需要更多的未知变量u i 。

所以未知变量个数随着基底函数阶数的增大而增大，基底函数的个数也随着阶数的增大而增大。

对于拉格朗日线性基底函数，在任何单元中都是由两个基底函数组成的：
1φξφξ==- (16)
这里ξ是基元中的本地坐标。

所以对于N 个单元，就有2N 个基底函数i φ。

二次拉格朗日单元有三个基底函数：
(1)(12)4(1)(21)φξξφξξφξξ=--=-=- (17)
所以对于拉格朗日二次单元，N 个单元有3N 个基底函数。

回顾前面讲的，方程(14)必须对所有的v ∈V 满足，下面我们考虑所有由基底函数线性组合构成的函数空间，即：
()()i i i v x v x φ=∑
(18)
但是由于v 线性加入(14)，这就说明，如果以任意基底函数i φ作为v ，方程(14)
都成立，那么对于所有基底函数的线性组合(18)也成立，即所有v ∈V 都成立。

所以(18)等价于一个(k +1)N 个方程（每个i φ都有方程(14)）、(k +1)N 个未知变量（u i ）的系统，这里k 是基元的阶数（k ＝1线性，k ＝2二阶）。

下面解释为什么COMSOL Multiphysics 有限元方法具有这样的作用。

COMSOL Multiphysics 建立(k +1)N 个方程(14)的集合。

首先，注意到Г(u )和F (u )都是u 的通用非线性函数。

所以一般没有闭合解。

在第一章中，我们介绍了COMSOL Multiphysics 中预置的零维问题非线性求解器，即f (u )=0，u 是唯一未知变量。

非线性求解器基于牛顿方法。

牛顿方法的N 维矢量方程为：
()0L U = (19)
这里U 是未知数u i 的矢量，L (U )是取代(14)中v 的基底函数i φ的方程系统，即：
000()()()K U U U L U -= (20)
这里的K (U 0)是刚度矩阵，L (U 0)是荷载向量。

刚度矩阵是L 的负雅克比矩阵：
00()()L K U U U
∂=-∂ (21) 所以(20)是给定近似解U 0情况下U 的线性方程。

因此，如果U 0能够足够接近真实解，线性方程(20)会找出更加精确的解U ，该过程可以不断重复，直到最后解足够精确。

显然，牛顿法非线性求解器是COMSOL Multiphysics 偏微分求解器的核心。

COMSOL Multiphyiscs 可以实现偏微分方程有限元分析的所有步骤。

如果可以的话，它用符号组成非线性运算符L (U )的雅克比矩阵，如果不行的话，采用数值方法建立雅克比矩阵。

如果偏微分方程本身是线性的，就不会很复杂。

但是建立刚度矩阵是一个非常庞大的工作——这在有限元分析中很常见。

不久以前，有限元分析，包括网格划分和建立刚度矩阵还经常是很多博士的研究课题。

对于新的偏微分方程组合，甚至变系数问题（例如1.2节介绍的准线性系统），建立刚度矩阵的工作量大的吓人。

进一步讲，方程(20)的稀疏求解方法和时间积分需要通过编程工作来调整使其适合于某一问题。

好在COMSOL Multiphysics 将这个工作作为一个子程序（MATLAB 函数）来完成，实现与复杂偏微分方程系统的无缝衔接。

在介绍边界条件实现方法之前，我们先通过一个常微分方程练习来强调一下刚才讨论的概念。

通过该练习了解一下弱形式，找出给定二阶常微分方程的近似数值解。

练习2.5
有限元详细计算实例。

使用有限元方法核心变分原理求解一个简单的常微分方程，详细介绍计算过程中出现的概念。

问题比较简单，二阶常微分方程如下：
22248d u u x dx += (22)
边界条件为：
(0)()04
u u π
== 使用了弱方程式。

该二阶常微分方程解析解如下（试证明之！）：
22()cos(2)1sin(2)218u x x x x π⎛⎫=+-+- ⎪⎝⎭ (23)
因为已经知道了解析解，就可以比较近似解的误差。

建立弱方程式的第一步是假设一个权重函数和试探函数。

将U (x )作为试探函数，()x φ作为权重函数。

我们一会儿再讨论这些函数的准确形式。

试探函数U (x )构成了常微分方程的解。

所以如果将U (x )带入(22)式，可以得到残差为：
22248d U R U x dx =+- (24)
接下来的工作就是最小化残差。

要最小化残差首先要计算权重残差，将方程(24)两侧同乘以()x φ，然后在整个域中进行积分（即0≤x ≤π/4）
2/4
220()48d U R x U x dx dx πφφφ⎛⎫=+- ⎪⎝⎭⎰ (25)
使用分步积分法简化上式为：
/4
/4200()48d dU R x U dx x dx dx dx ππφφφ⎛⎫=-- ⎪⎝⎭⎰⎰ (26)
为后续计算方便，我们需要做一些关键假设。

因为只要满足边界条件，我们可以任意设定U (x )和()x φ，我们令U (x )＝()x φ。

这被称作Galerkin 方法。

如果U (x )≠()x φ，给出Rayleigh-Ritz 公式。

我们需要选择一个x 的代数函数以满足边界条件u (0)=u (π/4)=0。

11220()()()N N N i i i x U x x u u u u ϕφϕϕϕϕ====+++=∑L
(27)
假设N 个方程为：
212444N N x x x x x x πππϕϕϕ⎛⎫⎛⎫⎛⎫=-=-=-
⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭L (28)
所以 1()4N i i i x u x x πϕ=⎛⎫=- ⎪⎝⎭∑ (29)
这种选择满足边界条件，且不需要考虑方程包含的项数。

由于U (x )＝()x φ，权重残差为：。