有限元方法的數學理論-變分問題簡介

04-06

有限元方法是由R.Courant於1943年首次提出並在20世紀50年代由航空結構工程師所發展起來的，我國馮康教授對有限元的理論基礎做出了巨大貢獻。有限元方法經歷了70年的發展已經是一個相對成熟的領域了，但在鞏固其物理數學基礎方面，以及擴大其應用領域，例如非線性、多體機構分析、材料微觀結構計算等，有限元方法依然會不斷發展。其實有限元線性剖分、離散化未知函數的想法是很符合我們直覺的，但是在嚴格化它理論基礎的時候依然是需要謹慎對待的，因此本文參考了《有限元方法的數學理論》一書，並加入自己的一些理解，旨在整理和介紹有限元方法的嚴格化理論。

一、變分問題

考慮兩點邊值問題： $egin{equation} left{ egin{aligned} &-u(x)+sigma u=f(x),xin Omega=[0,l] \ & u(0)=0,u(l)=0 end{aligned} ight. end{equation}$

若u為上述微分方程的解，對於一個滿足v(0)=0的足夠階可微函數v，我們在函數內積空間上有

$<f,v>&:=int_{0}^{l} f(x)v(x)dx=int_{0}^{l} [-uv+sigma uv]dx\&=-uv|_0^l+int_{0}^{l} [uv+sigma uv]dx\&=int_{0}^{l} [uv+sigma uv]dx:=B(u,v)$

此時我們記 $V:={v|v,vin L^2[0,l],B(u,v)<+infty ,v(0)=0}$ ，其中 $L^2[0,l]$ 表示[0,l]上的2次勒貝格可積函數，其第一類範數定義為 $||v||_{1,Omega }:=left[ int_{0}^{l}(|v(x)|^2+|v(x)|^2)dx ight]^{frac{1}{2}}$ ，則兩點邊值問題可描述為：

求 $uin V$ ，使得 $B(u,v)=<f,v>,forall vin V$ .

此表述稱為兩點邊值問題的變分問題或弱形式。變分問題的解稱為兩點邊值問題的弱解，如果函數 $u(x)in C^2[0,l]cap C^1[0,l]$ ，則u為兩點邊值問題的古典解。弱解所在的空間稱為容許（試探）空間，由於變分問題必須對V中的任一元素v成立，故稱V為檢驗空間，若容許空間與檢驗空間重合，稱其為能量空間，事實上容許空間與檢驗空間未必重合。

此類變分問題被稱為Galerkin變分問題，下面是相應的Ritz變分問題：

作如下二次泛函： $J(v)=frac{1}{2}B(u,v)-<f,v>$ ， $J(v)$ 的極小值問題為：

求 $uin V$ ，使得 $J(u)=min_{vin V}J(v)$ 稱為兩點邊值問題的Ritz變分問題。

由此可見，相較於我們討論Galerkin變分問題在推導分部積分時應用了邊界條件，Ritz變分問題對於邊界條件 $u(l)=0$ 並沒有作要求，這也是兩種變分問題的主要區別。

設V是一個Hilbert空間，稱B( , )為V上的對稱、有界、V-橢圓性雙線性泛函，若 $exists M,alpha in (0,+infty )$ ，使得：

$B(u,v)=B(v,u);|B(u,v)|leq M||u||_V||v||_V;B(u,u)geq alpha ||u||_V^2,forall uin V$ .

定理1：設V是一個Hilbert空間，B( , )為 $V imes V$ 上的對稱、有界、V-橢圓性雙線性泛函，f是V上的線性連續泛函，兩類變分問題：

求 $uin V$ ，使得 $B(u,v)=<f,v>,forall vin V$ & 求 $uin V$ ，使得 $B(u,v)=<f,v>,forall vin V$

有結論：

(i)任一問題有解，則解不多於一個；

(ii)任一問題的解，必是另一問題的解；

(iii)如果 $u^*$ 是它們的解，有性質： $J(v)-J(u^*)=frac{1}{2}B(v-u^*,v-u^*),forall vin V$ .

此定理給出了Galerkin變分問題和Ritz變分問題的等價性。

二、Galerkin逼近

設 $V_hsubset V$ 為V的有限維子空間(或稱為有限元子空間)，則Galerkin變分問題的Galerkin逼近為

求 $u_hin V_h$ ，使得 $B(u,v)=<f,v>,forall vin V$ .

令 $dimV_h=M$ ，其基函數係為為 ${varphi _i(x)}_{i=1}^M$ ，有

$u_h=Sigma _{i=1}^MU_ivarphi _i(x)$ ； $v_h=Sigma _{i=1}^MV_ivarphi _i(x)$

此時我們將求解未知函數的問題轉化為求解已知函數系的係數的問題，將它帶入變分問題的Galerkin逼近：

$&B(Sigma U_ivarphi _i,Sigma V_jvarphi _j)=<f,Sigma V_jvarphi _j>\&Sigma _{i,j}U_iV_jB(varphi _i,varphi _j)=Sigma V_j<f,varphi _j>\& extbf{U} extbf{K} extbf{V}= extbf{V} extbf{F}\&forall extbf{V}in V_h, extbf{K} extbf{U}= extbf{F} holds$

故有限元方程組(整體方程)為：

$extbf{KU}= extbf{F}$ ，其中

$extbf{U}=(U_1,...,U_M)^T$ ， $extbf{K}=(k_{ij})_{M imes M}=B(varphi _i,varphi _j)=int_{0}^{l} [varphi _i(x)varphi _j(x)+sigma varphi _i(x)varphi _j(x)]dx$ $extbf{F}=(F_1,F_2,...,F_M)^T=<f,varphi _i>=int_{0}^{l}f(x)varphi _i(x)dx$

K為(總)剛度矩陣，F是載荷向量。

現在我們取未知的基函數係為一系列的分段線性插值函數，即將定義域用網格劃分來近似方程的解，此時我們將域[0,l]分為N個單元(非等間隔)，記 $I_i=[x_{i-1},x_i],1leq ileq N$ ， $x_i$ 為節點，並記 $x_0=0,x_N=l,h_i=x_i-x_{i-1}$ ，取基函數系

$egin{equation} varphi _0(x)=left{ egin{aligned} &frac{1}{h_1}(x_1-x),0leq xleq x_1 \ & 0,xin else end{aligned} ight. end{equation}$

$egin{equation} varphi _i(x)=left{ egin{aligned} &frac{1}{h_i}(x-x_{i-1}),xin I_i\&frac{1}{h_{i+1}}(x_{i+1}-x),xin I_{i+1} \& 0,xin else end{aligned} ight. end{equation}$

$egin{equation} varphi _0(x)=left{ egin{aligned} &frac{1}{h_N}(x-x_{N-1}),xin I_N\ & 0,xin else end{aligned} ight. end{equation}$

易知 $varphi _i(x_j)=delta _{ij}$ ，在此基函數系下，我們有

$K_{i,i-1}&=B(varphi _i,varphi _{i-1})=int_{0}^{l}[varphi _i(x)varphi _{i-1}(x)+sigma varphi _i(x)varphi _{i-1}(x)]dx \&=int_{x_{i-1}}^{x_i}left( frac{1}{h_i}cdotfrac{-1}{h_i}+sigma frac{x-x_{i-1}}{h_i} cdotfrac{x_i-x}{h_i} ight) dx=-frac{1}{h_i}+frac{1}{6}sigma h_i\&i=2,3,...,N$

同理： $K_{i,i}&=frac{1}{h_i}+frac{1}{h_{i+1}}+frac{1}{3}sigma( h_i+h_{i+1}),i=1,2,3,...,N-1$

$K_{i,i+1}=-frac{1}{h_{i+1}}+frac{1}{6}sigma h_{i+1},i=1,2,3,...,N-1$

$K_{N,N}=frac{1}{h_{N}}+frac{1}{3}sigma h_{N}$

$F_i&=<f,varphi _i>=int_{0}^{l} f(x)varphi _i(x)dx=int_{x_{i-1}}^{x_i}f(x)frac{x-x_{i-1}}{h_i}dx+ int_{x_{i}}^{x_{i+1}}f(x)frac{x_{i+1}-x}{h_{i+1}}dx\&=int_{x_{i-1}}^{x_i}(f(x_i)+O(h_i))frac{x-x_{i-1}}{h_i}dx+ int_{x_{i}}^{x_{i+1}}(f(x_i)+O(h_{i+1}))frac{x_{i+1}-x}{h_{i+1}}dx\&=frac{h_i+h_{i+1}}{2}f(x_i)+frac{1}{2}[h_iO(h_i)+h_{i+1}O(h_{i+1})],i=1,2,...,N-1$

$F_N&=frac{h_N}{2}[f(x_N+O(h_N)]$

由於每個基函數只與前後兩個基函數有重疊部分，總剛矩陣是三對角的，我們得到了結構如同 $left( egin{matrix} * &* \ * &* &* \&* &*&*&\&...&...&...&\&&&* &* end{matrix} ight) left( egin{matrix}&U_1 \&U_2\&U_3\&...\&U_Mend{matrix} ight) =left( egin{matrix}&F_1 \&F_2\&F_3\&...\&F_Mend{matrix} ight)$ 的方程，由 $K_{i,i+1}=K_{i+1,i}$ 可知，總剛矩陣是對稱的，並且當劃分是等間隔的時候， $K_{i,i+1}=K_{i,i-1}$ .事實上，我們這裡求得的係數就是待求試驗解在節點上的取值，而基函數系決定了我們假定它的分布是分段線性的。總剛矩陣可以不是如上分布的，它是三對角的當且僅當劃分是順序編號的，若編號不同，稀疏矩陣的非零元素的分布也是不同的。稀疏矩陣非零元素在矩陣中分布的最大寬度稱為帶寬，而在電腦中對於稀疏矩陣的存儲只存非零元素，因此帶寬與電腦計算的上限是直接相關的。事實上，Matlab中可以處理的最大矩陣取決於電腦的內存，Windows系統留給Matlab的連續存儲空間大約800MB，也就是說若不採取順序編號，電腦只能計算30000×30000的矩陣。

三、誤差分析

由之前的命題，對於Galerkin變分問題和其Galerkin逼近我們分別有：

$forall vin V,B(u,v)=<f,v>;forall v_hin V_hsubset V,B(u_h,v_h)=<f,v_h>$

因此對於 $forall v_hin V_hsubset V$ ，根據B的線性性，有 $B(u-u_h,v_h)=0$ .

另外，由柯西不等式我們有

$&B(u,v)=int_{0}^{l}[u(x)v(x)+sigma u(x)v(x)]dx \&leq max { 1,sigma}int_{0}^{l}[|u(x)|cdot|v(x)|+|u(x)|cdot|v(x)|]dx \&leq C_1left[ int_{0}^{l}[|u(x)|^2+|u(x)|^2]dx ight] ^{1/2}left[ int_{0}^{l}[|v(x)|^2+|v(x)|^2]dx ight] ^{1/2}\&=C_1||u||_{1,Omega }||v||_{1,Omega }$

$B(u,u)=int_{0}^{l}[(u(x))^2+sigma u^2(x)]dxgeq min{1,sigma } int_{0}^{l}(|u(x)|^2+|u(x)|^2)dx=C_2||u||_{1,Omega }^2$

結合上述性質，我們有

$&C_2||u-u_h||_{1,Omega }^2leq B(u-u_h,u-u_h)=B(u-u_h,u-v_h)leq C_1||u-u_h||_{1,Omega }||u-v_h||_{1,Omega }\&||u-u_h||_{1,Omega }leq frac{C_1}{C_2}||u-v_h||_{1,Omega }$

定理2：若 $u&u_h$ 是Galerkin變分問題和其Galerkin逼近的解，則下列結論成立：

(i) $B(u-u_h,v_h)=0,forall v_hin V_h$ ；

(ii) $||u-u_h||_{1,Omega }leq gamma inf_{v_hin V_h}||u-v_h||_{1,Omega },where gamma is a constant.$

可見，在範數 $||*||_{1,Omega }$ 下，可以將B(*,*)看做是向量相對於多重線性映射的內積，近似解 $u_h$ 為弱解在有限元子空間 $V_h$ 上的正投影，且近似解的誤差 $||u-u_h||_{1,Omega }$ 分析等價於函數插值逼近的誤差分析。

任一函數u的插值定義為 $ar{u}=Sigma_i u(x_i)varphi _i(x)=Sigma_i u_ivarphi _i(x)in V_h$ ，特別地，我們有：

$||u-u_h||_{1,Omega }leq gamma ||u-ar{u}||_{1,Omega }=gamma left[ Sigma _i||u-ar{u}||^2_{1,I_i} ight] ^{1/2}$

記： $e_i:=(u-ar{u})|_{I_i},|e_i|_{n,I_i}=left[ int_{x_{i-1}}^{x_i}|e^{(n)}_i(x)|^2dx ight]$ ，就有 $||e_i||_{1,I_i}^2=|e_i|^2_{1,I_i}+|e_i|^2_{0,I_i}$ .

因此接下來我們將要估計每個區間上的 $|e_i|^2_{1,I_i}&|e_i|^2_{0,I_i}$ .鑒於這裡的過程較為繁瑣(事實上需要藉助Sobolev空間嵌入定理確定誤差 $e_i(x)$ 的光滑性)，我們只給出結果，記 $h=max{h_i}$ ，有下述定理：

定理3：若u屬於(2,2)階Sobolev(Banach的)空間 $H^2(Omega ):={v|D^alpha vin L^2(Omega ),forall |alpha|leq 2}$ ，其中D表示廣義導數， $L^2(Omega )$ 表示2次勒貝格可積函數空間， $V_h$ 為分段線性插值函數系張成的有限元子空間， $u&u_h$ 是Galerkin變分問題和其Galerkin逼近的解，則

$||u-u_h||_{1,Omega }leq Ch|u|_{2,Omega },where C=gamma sqrt{1+h^2} .$
推薦閱讀：

TAG:有限元分析FEA | 計算數學 | 計算理論 |