一类带复合边界条件的偏微分方程的有界和有界区域的有界解.docxVIP

一类带复合边界条件的偏微分方程的有界和有界区域的有界解 考虑到带复合边界条件的偏微分方程。 {-αΔu+βu=s,x∈Ω,u=d,x∈D,α?u?ν+γu=r,x∈Γ,(1) 其中Δ为拉普拉斯算子,?u?ν为u沿单位外法线方向的方向导数,Ω为有界区域, 边界D,Γ满足D∪Γ=?Ω. 系数α、β、s、d、γ、r为给定常数, 在实际计算中, 这些系数可以是给定的函数. 这类方程在生物数学中常用来描述微环境因子(例如氧气)反应和扩散过程的拟稳态, 其中Δu表示扩散项,βu(连同常数项s)表示反应项. 本文以二维情形为例, 推导此类问题的有限元解法的代数方程组形式, 并给出一种快速计算其刚度矩阵的方法. 对于边界条件的处理请参阅文献. 1 单元区域中的真解 用有限元方法近似求解偏微分方程的数值解, 基本想法是把场域(例如问题(1)中的Ω)分割成很小的子区域,通常称为”单元”或者”有限元”, 然后对每个单元进行独立的处理和运算, 通过选取恰当的尝试函数(局部上与真实函数近似的简单函数), 使每个单元的计算变得简单, 经过对每个单元重复而简单的计算,再将其结果总和起来,便可以得到用整体矩阵表达的整个区域的解. 如图, 本文对场域采用三角剖分, 在每个单元(e)上, 真解可以用一个简单的平面方程ue来近似, 它是由3个节点处的尝试函数{φei,φej,φek}线性组合而成 ue=ciφei+cjφej+ckφek. (2) 为简单起见,以下省略尝试函数的上标e, 其选取如下(仍是平面方程) φi=αi+βix+γiy, (3) 使得 φi(xi,yi)=1,φi(xj,yj)=0,φi(xk,yk)=0, (4) 节点j、k处的尝试函数φj和φk的选取类似, 形式也完全形同. 将条件(4)带入函数(3), 通过解一个三阶线性方程组得到{αi,βi,γi} 的值,从而可确定节点i处的平面方程φi. 同理可确定节点j、k处的平面方程φj和φk. 2 单元n,l 接下来在弱解的意义下推导问题(1)的有限元解法的代数方程组形式. 假定u是问题(1)的解,则u在任意单元(e)上也满足此方程, 记该三角形单元为Ωe,对于任意在Ωe上具有紧支集的函数v∈H10(Ω)满足 ∫Ωe(-αΔu+βu)vdσ=∫Ωesvdσ, (5) 由格林公式, ∫Ωe(α?u·?v+βuv)dσ=∫Ωesvdσ, (6) 令u=3∑i=1ciφi,v=3∑j=1bjφj,带入式(6)整理得 3∑j=13∑i=1∫Ωe[α(?φi?x?φj?x+?φi?y?φj?y)+βφiφj]dσcibj=3∑j=1∫Ωesφjdσbj.(7) 由v的任意性,即bj的任意性知 3∑i=1∫Ωe[α(?φi?x?φj?x+?φi?y?φj?y)+βφiφj]dσci=∫Ωesφjdσ(8) 对j=1, 2, 3 成立. 写成矩阵的形式即 (αA+βB)c=sb, (9) 其中向量c=(ci)3×1为待求量,A=(Aij)3×3,B=(Bij)3×3,b=(bj)3×1各分量表示为 Aij=∫Ωe(?φi?x?φj?x+?φi?y?φj?y)dσ,Bij=∫Ωeφiφjdσ,bj=∫Ωeφjdσ.(10) 由前所述, 若对每一个单元(e), 均已求出了所有如式(10)所示的各分量, 再将各单元的局部矩阵组装成整体的矩阵(刚度矩阵), 并联系边界条件删去相关节点, 得到最终的线性代数方程组, 就可以根据矩阵的规模选用合适的方法求解方程组了. 所以快速生成刚度矩阵就成为关键. 3 bij和bj的计算方法 定理1 式(10)中的分量可以按以下方式计算, Aij=(βiβj+γiγj)S,Bij={S/6,i=j,S/12,i≠j,bj=S/3,(11) 其中:i,j=1, 2, 3;S是三角形单元(e)的面积. 证明由于选用平面方程作为尝试函数, 其表达式具有式(3)所示的简单形式, 计算Aij是容易的.下面给出Bij和bj的计算方法. 不失一般性,设单元(e)为等腰直角三角形区域,边长为1,从而面积S=1/2.以节点i为原点建立如图2所示的坐标系, 平面φi如图2中阴影所示,φj,φk与其类似,其平面方程依次为 {φi=1-x,φj=x-y,φk=y. 直接计算可得 ∫Ωeφidσ=∫Ωeφjdσ=∫Ωeφkdσ=16=S3. 同样∫Ωeφ2idσ=∫Ωeφ2jdσ=∫Ωeφ2kdσ=112=S6, 以及∫Ωeφiφjdσ=∫Ωeφiφkdσ=∫Ωeφjφkdσ=124=S12. 注意到式(10), 从而得到式(11). 在计算实践中, 按照式(11)计算刚度矩阵要比直接计算式(10)中的重积分或其他数值积分方法如蒙特卡罗法明显要快.