偏微分方程的有限元法求解.docx

偏微分方程的有限元法求解.docx

《偏微分方程的有限元法求解.docx》由会员分享，可在线阅读，更多相关《偏微分方程的有限元法求解.docx（17页珍藏版）》请在冰豆网上搜索。

偏微分方程的有限元法求解

16.901讲义笔记

一维有限％

首先，我们考虑•个比上一节稍微复杂点的问题;豎二f（X）,卫冲，V（O）=O.V（L）=O

在这里，f（X）是）C的般函数，我们来看•个特别的情形：

f（x）=x（L-x）,此时，方程的梏确解如F：

有限元方法利用加权残差的方法■其中:

（1）

设va）=£«Ma）,v（）（）是我们对v（x）的近似，

省为未知常数9V|（x）是用户选择的歯数，即形状朗

数：

（2）

定义N个加权残差

L

Rj=p^（x）R（V）dx•j=l->Nt

o

其中，RV）二器・f为绒差

凹⑴足“用户”选择的加权函数，即权函数:

（3）

令加权残并为冬•町以确定⑷的值，即求耳使得对所fi1=I->N,Rj=Oe

令限元方法（）是加权残若法的一种，下血看看

我们是如何用它来解决问题的。

一维有限元方法

有限元方法（〉扌野个连续区域离散化-系列小

单尤，这些单元与有限差分法（）或有限体积法

（）产牛的网格完全相同，而佼之前两者主耍的优点在于：

能够容易地把握单元的变化范囤。

对于我们讨论的一维问题,可以将区域（数轴〉离散化为如下图所示:

这里，叫三单•元的个数。

我们还会用別下血i些定义：

个三角划分；尽管令限元法对于一维，二维，三维甚至高细问题都是仃效的，们我们还是要谈及区域离散化的一种方浓，即三角划分。

4T定义为第I个单元所在的区域。

对于_维问题,

这表明，TS-个满足片心的X的集合。

接卜来耍确定的是毎个单兀该用什么样的函数，典型的函数形式就是用从一个单元到卜一个单兀保持解连

续的多项式。

例如：

一个线性有限元如卜團;i示:

在毎个单元内的函数是线形的，在毎两个单元的交点处足连续的。

对于专门诜择的满足线件变化的形状函数,右估计残差时有一个很明显的问题：

回忆前曲的内容，RV）二器一f,它在一个单冗里等于什么呢?

因为函数是线性的,所以器=0,则有:

R（V）=f,即R（V）与无关。

冋时，满足线性变化的形状函数似乎也是一个好的近似，我们举-个例子来说明。

如果我们把-个区域三角划分为10个相同的小区域，采舟线性插值求苴精确解，这是一个合理的解。

（见附图〉

实际上，的标灌方法足利用分部积分浓〔），

将対V求导转化为刈抑求导：

Rj=

在这种形式下，我们能用线性的形状函数，这就足标准的方法。

（如杲可能的话，分部积分法对V和

具有相同阶数的导数。

）

注意，*分部积分汉的另一个优点足保讦遵循所观测到的导怛律（尽骨这个优点不是很容易农现）

*在积分和微分时我们嗖格外小心•因为在节点处，分段线性函数的导数不存在，这也足一个高级的话题。

节点基

存很多种不冋的方式描述网格」；的分段线形凶数，在有限元法屮，通常使用的是所谓的“节点基”。

V（X）的节点基的形式如下所示：

V（X）=£Wg（X）,这里.％是未知常数，

§（X）娈满足以下条件：

尙（卸=°"j在同一个节血处

所以，对线性单元,我们有;

由此,我们注意到:

V（斗V衙（申_片英中，5（片）除了i=j时都为

所以,片是一个常数,表示节点处的解。

•绒参考元和内无式|

在一个单尤屮，我们需耍为两个节点构适线性内井

K+歹而三在单元丁内节点i的节点插值，

G+新在甲元壬内节点i+1的节点插隹尽■维线性元内节点插值式可容易地表为:

但対十高阶的（例如二次，三次等）、尤瓦是直绒的单元,这是非常麻烦的。

所以,我们耍做的是把单元映射为主单元或参考单元。

-维参考元

将x轴上以X|和甩为节点的单元映射刘§轴上节点为&严+1和©=+1的参考单元上。

在参考单元里，线性内插式:

为：

映射过程保持歹一>X,反z亦然。

考虑到多纽的单元，我们将用等参数映射。

特别地，我们设对给定的：

位置，x位置可表为：

x（5）=粘忆）+妬点）或者刈計=右船）显然，有；

X<,）=X-1）=X因为<（T）=1，G（-i>=0x（：

）=Ki）=X：

因为以⑴i,5

（1）=0

在中点，^=0»我们有：

X0）=&（o）十Xc,

2（o）=+

x（0）=卅乜）为平均值

所以,一般函数的内差式（例如函数对应方程的解）在给定的•维线性元内拥冇下血的形风；

/需耍求出W

▼⑴帚专砖似X））

对于-•维线性元,我们能很容易地衍到逆映射;

-幣理解得c

5（x）=~r2

扣：

7）

注意:

容易验证眞片）=-1利§（均）=1。

所以.为了右三角单元中找到v（x）,我们可以采取以下步骤：

（1）利用員狗映射找到g的位宣

（2）计算形状苗数在占的值，例如:

<3）估计va）=隔（$）+«2⑷

利用参考兀计养梯度

为了进行有限尤分析，我们通當需耍在一个单元中仙

计v（x）的梯度，利用参島兀的等参数映対冇:

密卜仝V|f|G（切+%2（乱談））

再用他式法则，可展开为：

毛煌卷

回忆前面的维参考元;

d£]

x-£d筑x）二

乂曲扣f'

将上面几项带入原等式，有:

dV

（不兄在边界上）

回顾有限元方法的加权残差衷迖式:

除了在边界*=0和x=L1:

7只有第二和第三项可能是非零项。

利用伽辽金冇顒兀汕计舁i兽兽企

伽辽金加权殁差広即鬥=®，所以〜不为零，其值在节点j的周帽单元内变化，例如：

TH.Tr

200心

竺1岂dxIdxdx

；d®dV.『d®dVj

珂F詰".*不如现在.我们可以展开：

v（x）龙％a匸讐龙v.警

0I10

所以，例如：

cbc

Aidxdx台£*

但现在在三角单元0内,只有上j-hi=j对应的岂糾不为零！

dx

［吆岂心训丄吆虹却叭吆业$

*||dxdx1知dxdxFidxdx最后，因为&Jj»则有:

1吆岂关畛如％生1业改

•"ndxch1-

剰余的步骤就是利用我们的丰单元内斧式和映射找岀

每一项二

瞥2焙訣

这时利用链式法则，可得；

「芝j寸［埜岂［虫生［dx

J*dxdx丄l码dx哲dxJ

这是把对x的枳分化为对£的枳分，则上式变为:

!

*此、％4_制卍；型\（丄生血1盂詁叮〔苕訓灵亚①这里，我们已经用过d“会改或者曲T±l蛮和我们前面提到的一样：

°»1・$产+1

d<.I

17=~21?

=2东_岭_齐

nr審訣訂店占界占）（宁卜=-2（^）r>

一几一1

所以，用类似的方法我们可以得flJK他的非咨项:

而口，对丁线性元,我们有:

计算禹fdx

0

丙为fg可能是X的一个普通歯数，我们阳亥用一般的积分技巧來求bjfdx.至少可以足很好的近似。

0

和前面一样.除了在单元W内，叫=0

L

=>（網fdx=[少fcbc4[Q]fdx

我们还像从前一样利用参考兀将対X的积分转化为

对石的积分,

仲借）吗[

无论哪种情况,

我们都会碰到一个潜衣的困雄，就是

卜面这种形式的一维积分M

Jfl（^）df这里的g（即是£的普通函数

I高斯求柄

高斯求积就是把一个积分近似列被积函数曲数值的加权求和，即；

厂IN

J,g⑷如工叩⑷

这里的厲.疳是通过对多项式g（0耕确积分来确定的，看下而这个具体的例子，

刈于所有q的值，OQSM

首先,我们注意到:

1-1I』>J

N“：

_亦斯求积法则

我们遇过单个加权点寻找最高阶多项式的积分;f,g&）d“qg（G

2q十*q十*G.十…=⑷{$十十鸟济'十G铲十q昇十…}

现在，逐顼进行比较有：

G：

2=qn|“产2|

q:

O=a^=2^,=>（£=0|

G.:

^a^z=0=>不满足

=>

j'g（^^2g（0）

对于线性多项式,能楮确求积

=2髙斯求积法则

接下来,我们用两个点

2%吗屮典"“=剣勺吟Mi略喇引吗何斤％彳,

g（“I即+口2冷）+C4（°l匸+al^z”…

逐项对比有:

q;2=q十冬

q：

O=a&+a屋

Q：

0=a1^+o,^-

2

G:

y=^i+a2^2

对这些约束条件，我们可以采取如卜操作：

a丘=-a.fi

a^3=-a^x

»厂垛假设关于是对称的

此时，粥珅q有：

-4

再看对于昭的约束：

«!

+tr:

=2

2购=旳=>!

«,=g,=l\

最后，来看对于G的约束:

g（f）df=g（^^-）4g（卜书）

这个公式对丁-三次多项代也是精确的，例如:

的求积过程中，我们发现:

其中g（G・W）f@）（¥）

边界条件的殛

方程^=f（x）（或Vv=f）的边界条件通常有卜価的两种；

1，狄里希茱条件：

V在边界上有确定的值；

2,纽曼条件:

里在边界上有确定的值

dx

处理这些边界条件的标准方法如下：

狄里希莱条件；

♦要使V在边界上等于一个给定的值,并且令所有的权甫函数在边界上为零，则加权残差方程幻

在H前我们止在研究的问题中，如果令V/|-x=O,x=L时为零,有：

对所有的j,"j（0）=0,®（L）=0

RV（0）=0.V（L）=O

注盘：

V（0）=£vq（O）"=O除了i=l.RO）均为零

II

V（L）=Jvi^L}=V^,=O

II

如果我们比较方程和未知数的个数，会白：

未知数：

V,i=l->Ne+l

R）=0|=2-*Ne

方程：

v,=o

V^.I=0

=>N,+I个未知数

o

•=>N.4-I个方程

纽址条件:

「假使令款Z时为》即孰）小

*最后,除了y=i时，叫（o）=i,恒有^（0）-0；

Wilt,>1对应的加权残苏的形式为；

*加权妓丼边界约束具有如卜形式：

RUf罢罢dx.V』罢罢妇鳥冲

和前面一样，一哄anc+i个未知数，而对血的方程:

也足叫J个，所以可以求解中

dWdxdM2dx

%对所冇单元循坏汁畀刚性和妙差

%

forii=l:

Ne,

knl=ii：

kri2=ii-H;

xl=x（knl）:

x2=x（kn2）;

dx-x2-xl:

dxidx=2/dx;dxdxi=l/dxidx:

dMldxi=-1/2:

（iM2dxi=1/2;

■dNldxi*dxidx;

=dM2dxi»dxidx;

%用高撕积分估计力项的积分tornn=1:

NGf,

、得到高斯点的xi

xiG=xiGH（nn）;

%求N1利N2（即在xiG的权重/插值〉

N1=0.5*

N2=0.5*

%对高斯点求f

fG=xG*（l-xG）：

%在节点处估计权凶数在斋斯点的被积负数gCl■Nl*iG*dx

rC2=N2*tG*dx

%利用上面结果来修正右端项b（kul）・b（knl）+aGf*（nri}*gGl;b（kn2）=b

end

end

%在x=0设置Dirichlet条件kn