应用PDE讲义15_RitzGalerkin法.pdf

1 应用偏微分方程与科学计算 讲义 十五 Lecture Notes on Applied Partial Differential Equations and Scientific Computing No 15 马 石 庄 2011 11 03 北京 2 第第讲讲Ritz Galerkin 方法Ritz Galerkin 方法 教学目的 教学目的 Ritz Galerkin 方法是解数学物理方程定解问题的经典方 法 由于逼近空间不同而与差分法截然不同 有限元法从变分原理出 发 把 Ritz Galerkin 方法与差分方法的有机结合 将逼近空间分割 成许多有限单元进行分片插值 从而具有高度的灵活性和广泛的适用 性 主要内容 主要内容 1 Ritz Galerkin 方法 3 1 1 Ritz 方法 1908 5 1 2 Galerkin 方法 1915 8 1 3 权余量方法 10 2 一维有限元 13 2 1 基于 Ritz 观点 16 2 2 基于 Galerkin 观点 19 2 3 逐单元计算 21 3 有限元法扩充 25 3 1 高次插值 25 3 3 混合有限元 27 3 1 基于变分原理的差分格式 28 练习 15 28 3 由于计算机运算的高速度 在古老的差分法和 Ritz Galerkin 变分 方法获得新的生命力的同时 另一种新的强有力的离散化方法 有限 元方法被创造出来 有限元法是经典的 Ritz Galerkin 变分方法与差分 方法的有机结合 传统的变分方法由于采用整体逼近空间而与差分法 截然不同 有限元法从变分原理出发 却将逼近空间分割成许多有限 的单元进行分片插值 从而具有高度的灵活性和广泛的适用性 有限元思想最早可以追溯到 R Courant 在 1943年的论文 1956年美国工程师M J Tuner R W Clough 等从结构力学角度重新提出后 在 20 世纪 60 年代逐渐形成系统的方法 1965 年 中国科学家冯康在 应用数学和计算数学 杂 志上发表题为 基于变分原理的差分格式 的论文 解决有限限元法 的收敛性 稳定件和误差估计等基本理论问题 从而在世界上首先了 有限元法完整的理论基础 这是计算方法发展史上的一个里程碑 也 是中国科学家对现代数学为数有限的贡献之一 1 Ritz Galerkin 方法 1 Ritz Galerkin 方法 先考虑两点边值问题 d d d d 0 0 其中 是一次连续可微函数 min 4 0 上的连续函数 0 平方可积函数 这里 平方可积函数空间 Soblev 空间 进一步定义 满足 0 满足 0 0 显然 构造二次泛函 作用积分 d d d d 和 1 2 建立与定解问题 在一定条件下等价的变分问题变分问题 求 使得 min 和变分形式变分形式 求 满足 在力学上 变分问题 相应于最小位能原理 变分形式 相 应于虚功原理 左端项表示虚内功 右端项表示虚外功 在定解问题 中 要求解 而在变分问题 和变分形式 5 中 解只需要具有一阶广义导数 即 且满足边界条 件 因为 是无穷维空间 直接计算是困难的 假设在其中找到一 个维数 的有限子空间 设空间 一组基函数 即 span 对于 中任意一个函数 都有 其中 都是常数 称 为试探空间 1 11 1 Ritz 方法 1908 Ritz 方法 1908 瑞士物理学家 Walther Ritz 1878 1909 使用变分问题 用试探函数逼近的近似问题是 求 使得 min 根据二次泛函 和 的性质 有 1 2 1 2 1 2 6 其中 是已知的基函数 所以 和 是 已知的实数 求 极小的问题就转化为求以系数 为变 量的二次函数的极值问题 在前面已经证明 在 维欧氏空间 引进向量和矩阵记号 表示括号内向量或矩阵的转置 令 定 义 的内积 考虑个变量的二次函数 它在 取极值的必要条件是 0 0 1 2 假定 即 为对称矩阵 则 7 2 不难看出 若令 1 2 则二次函数 关于 取极值的必要条件是 是线性代数方程组 的解 根据这个结果 变分问题 的解的充分条件为 0 0 1 2 即 0 0 方程有唯一解 值得注意的是 如果记 线性代数方程组 表示为 系数矩阵 8 一般是满秩的 需要进行计算大量的积分运算 此前 Rayleigh 以 类似的思路求解特征值问题 学界也称 Rayleigh Ritz 方法 1 2 Galerkin 方法 1915 1 2 Galerkin 方法 1915 1905 年 曾在中国参加中东铁路修建的俄国工程师 Boris Galerkin 1871 1945 使用变分问题 用试探函数逼近的近 似问题 求 满足 0 设 的任意元素 且变分问题的解为 代入 有 0 9 此式对任意的 成立 即对任意的 维向量 都成立 所以有 0 0 在 是对称的条件下 Galerkin 的结果 与 Ritz 的结果 相 同 如果 非对称 仍然可以求 Galerkin 变分问题 因 此 Galerkin 法可以看成是 Ritz 法的推广 应用范围也确实比 Ritz 法更宽 经常把两者合称为 Ritz Galerkin 方法 再说明以下几点 1 当边界条件非齐次时 由于非齐次的自然边界条件只影响 项 因此不影响 即 的选取 而对非齐次的强制边界 条件 如 则应取 这里 属于 且满足 的任一固定函数 然后再按照 前面的过程同样处理 得出类似的线性方程组 2 Ritz Galerkin 方法的困难在于坐标函数的选取 对两点边 值问题来说 坐标函数的选取困难比较容易 而对高维椭圆型方程的 边值问题来说 除了在一些非常规则的区域以外 选取满足强制性边 界条件的函数就是非常困难的 这是 Ritz Galerkin 方法得不到有效 推广的最大障碍 3 在用 Ritz Galerkin 方法导出线代数方程组 或 时 要 通过计算大量的积分来得到 和 项 一般说来 当 10 近似空间取为 维时 要计算约 次 非对称时 或 2次 对称 时 其计算量是相当可观的 4 即使形成了 阶线性方程组 在按传统方式选取坐标函数 时 方程组的条件数很大 计算时数值不稳定 另外内于系数矩阵是 满秩的 求解时所花费的时间和占用空间都是惊人的 5 Ritz Galerkin 方法的推导线性方程组的过程 对于将函数 集合取为所述的 是完全相同的 6 这个方法得到近似解的形式与有限差分法有一个明显不同 它得到的是精确解在所考虑的区域上的一个近似函数 或者说是一个 连续形式的解 而不仅仅是精确解在有限个结点处的近似值 从上述的第 2 4 点看到 有必要寻找一种好的方法 尤 其是选取坐标函数的方法 以保持 Ritz Galerkin 方法的基本思想而 弥补其缺陷之处 第 5 点又使得寻找这样一种方法成为可能 因 为以 中的函数的选取将比 具有更大的灵活性 1 3 权余量方法 1 3 权余量方法 有许多种途径可以导出有限元方法和谱方法 带权残量法 Method of Weighted Residuals 是其中之一 它的概念简单而数学 上严格 权余量法的名称是 Crandall 1956 引入的类似的想法 Collatz 1960 也考虑了用 误差分布原理 error distribution principle 的名称 Harrington 1968 在标题 矩法 method of moments 也有类似的概念 11 对于偏微分方程的边值问题 0 0 其中 是偏微分算子 求该定解问题的近似解 设近似解形如 其中 是已知的函数 称为试探函数试探函数 trial function 用 于满足初始和边界条件 系数 待定 令残量为 为了获得 的方程 令带权 残量的内积为零 0 1 正是以这样的方式求近似解的 都可以冠以带权余量方法 method of weighted residuals 权函数 称为检验函数检验函数 test function 显然 为了保证上述 个关系式独立 要求 是相互独立 在 0 1 意义下正交 如果 时 是一个完备函数系 则残量 一 定与 的每一个函数正交 这意味着 时 残量 按内积诱 导的范数 收敛于零 如果近似解满足初边值条件 那么可以期 望近似解 按范数 收敛到真解 即 lim 0 12 公式 的重要性在于把许多看来没有联系的近似方法联系起来 差 别在于选择不同的检验函数 如果选择检验函数就是试探函数 1 那么就是 Galerkin 方法的原型 1915 求 满足 0 Kantorovich 和 Krylov 1958 建立了 Galerkin 方法与 Rayleigh Ritz 方法的联系 转变为求变分 求 满足 0 配置方法 在区域 内规定的 个配置点 上 取检验 函数 为 Dirac 函数 则残量为 0 配置法的名字和第一个例子由 Frazer 等 1937 提出 1941 年 Bickley 应用配置法到抛物型方程 最小二乘法起源于 Gauss 问题是求 使得余量 在均方 意义下最小 min 直接可以取检验函数为 1 2 因此 13 0 等价与选择 使得 极小 2 一维有限元 2 一维有限元 除了在一些规则的区域外 要应用经典的 Ritz Galerkin 方法 选 择取满足边界条件的基函数常常很难办到 不但如此 系数矩阵一般 不是稀疏的 当方程组阶数较高时 矩阵元素的形成需要计算大量的 积分 计算量及的存储量都相当大 20 世纪 60 年代中期以来迅速发 展的有限元 提供了选取基函数的新思路 克服了传统的 Ritz Galerkin 方法的缺点 有限元法与传统的 Ritz Galerkin 法的 主要区别在于 应用样条函数方法提供了一种选取 局部基函数 或 分片多项式空间 的新技巧 从而在很大程度上克服了 Ritz Galergin 方法选取基函数的固有困难 有限元法实质上是 Ritz Galerkin 法的一种 有限元法的基本问题可归纳为 把问题转化成变分形式 选定单 元的形状 对求解域作剖分 构造基函数或单元形状函数 形成有限 元方程 提供有限元方程的有效解决 收敛性及误差估计等 只讨论一维情形 对模式问题 及其等价的变分问题 和 将给定的空间区域 剖分 设节点 满 14 足 得到子区间 0 称为 单元 单元的长度 且记 max 考虑区域 上的分段线性函数 在区域上连续 在每个子单元 上是线性函数 这些函数的全体构成集合 在每个 上是线性函数 这样定义的函数集合与剖分有关 因此有下标 标记 显然当 时 有 0 单元形状函数 由位置函数在节点上的一组值 1 1 15 而言 按照线性插值公式 试探函数为 0 为使计算过程标准化 引入仿射变换 0 把单元 变为 轴上的参考单元 0 1 在参考单元上 插值表达式简洁 定义两个线性函数 1 0 1 则线性插值公式写为 0 1 1 1 1 16 2 1 基于 Ritz 观点 2 1 基于 Ritz 观点 从 Ritz 观点出发 变分问题 近似问题是 求 使得 min 因此有 1 2 而 1 2 1 2 d d d 1 2 d d d 其中试探函数为 0 和 1 17 d d 0 则 d d d d 且 d d 其中 都不是积分变量 可以写为 的二次函数 1 2 2 其中 d d d 18 d d 合并同类项后 注意到 0 得到 1 2 其中 其中 定义 0 0 二次泛函 取极值的必要条件是 0 1 2 也就是关于 的三对角线性方程组 2 1 这里用一个性质 2 1 从而近似解为 T 满足 19 min 2 2 基于 Galerkin 观点 2 2 基于 Galerkin 观点 从 Galerkin 观点出发 变分问题 的近似问题是 求 满足 0 自然对每一节点构造山形函数 1 1 1 1 0 其它 而 1 0 其它 容易验证 1 1 1 20 因此 问题是 求 相当于求 使得 1 借助单元 上的仿射变换 及 轴上 0 1 上的标准山形函数 1 0 1 则对 0 其它 1 和 0 其它 显然 当 2时 0 以 为基函数 近似问题 表述为 21 0 1 注意到系数矩阵的第 行只有三个非零元素 即 d d d d d 对比得知 类似地 可以证明 和 再次见到三对角稀疏线性代数方程组 因此 无论从 Galerkin 观 点和从 Ritz 观点出发 都得到形式相同的有限元方程 区别在于 从 Galerkin 观点出发 建立有限元方程更直接 且更具普遍性 当 非对称时 无法用 Ritz 观点求解边值问题 2 3 逐单元计算 2 3 逐单元计算 在实际工程计算中 并非直接形成有限元方程 而是逐个单元分 析局部二次型即单元矩阵 称为单元刚度矩阵 element stiffness 22 matrix 再由单元刚度矩阵形成总矩阵 称为总刚度矩阵 global stiffness matrix 对于 Ritz 观点而言 对第 个单元 记 和 则有 把单元刚度矩阵 扩充为 阶方阵 0 0 0 0 在第 1行和第 行与第 1列和第 列交叉位置的元素为 的四 个元素 其余元素都是 0 记 T 则 T T 于是 23 其中 对于 Galerkin 观点而言 不是把方程组系数矩阵的每个元素 按 和 的顺序依次计算 而是按单元逐个计算 在累加起 来得到方程组的系数矩阵和非齐次向量 例 例 考虑边值问题 用线性元求解 d d 2 0 1 0 0 1 0 对应的 Galerkin 意义变分问题为 求 满足 d d d d d 2 d 0 将区间 0 1 分别划分为 4 个子区间 节点为为 0 1 4 1 2 3 4 和 1 从左至右依次编号为 显然 1 4 1 4 标准山形函数 1 0 1 求得系数 1 d 1 1 d 1 24 1 d 1 2 d 2 d 单元刚度矩阵为 1 1 1 1 1 0 4 单元荷载向量 1 1 0 4 依次扩大 1 1 1 1 1 0 0 0 1 0 1 1 1 1 0 0 1 0 0 1 1 1 1 0 1 0 0 0 1 1 1 1 和 1 1 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0 0 1 1 得到 25 1 1 1 11 1 1 11 1 1 11 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 和 1 1 1 1 1 1 1 1 1 2 2 2 1 考虑到约束条件 0 在总刚度矩阵 其中划去第一行和第一列 相应地 在 中也划去第一行 得到有限元方程相应的代数方程组 4 2 1 12 1 12 1 11 1 4 2 2 2 1 求得解为 7 16 12 16 15 16 1 3 有限元法扩充 3 有限元法扩充 3 1 高次插值 3 1 高次插值 以上的讨论都是在单元上进行线性插值的 得到分片线性函数的 26 试探函数空间 如果希望提高逼近精度 一个有效途径就是单元上的 形状函数取为高次多项式的情况 要在单元上确定一个 次多项式 应 该给出 1个条件 除了前面叙述过的 在 和 上的函数值 和 之外 还要给出 1个条件 有下述两种主要的给法 将给定的空间区域 剖分 设节点 满足 在子区间 0 上讨论插值问题 单元的长度 为使计算过程标准化 引入仿射变换 1 2 0 把单元 变为 轴上的参考单元 1 1 称 为局部 坐标 1 1 为标准单元