《計(jì)算方法》課件2第5章

5.1邊值問題的變分形式5.2瑞茲－伽略金方法5.3有限元法

5.4實(shí)例介紹習(xí)題五

5.1.1二次泛函的極值

首先，我們引入向量和矩陣的記號(hào)

5.1邊值問題的變分形式x=(x1,x2,…，xn)Tb=(b1,b2,…，bn)TA=

考慮n個(gè)變量的二次函數(shù)

在x0=()T取極值的必要條件是

(Ax,x)-(b,x)即

(aik+aki)-bk=0,其中k=1,2,…,n

如果A是對(duì)稱矩陣，即aik=aki，則上式化為

aki-bk=0,k=1,2,…,n

若令

(5.1)則J(x)在x0取極值的必要條件是x0是線性代數(shù)方程組Ax=b的解。

下面考慮實(shí)變量λ的二次函數(shù)

因?yàn)锳是對(duì)稱矩陣，則

(5.2)因?yàn)镴′(x0)=0，即對(duì)于任意x

φ′(0)=(Ax0-b,x)=0

即Ax0-b=0，x0是Ax=b的解。因?yàn)镴(x)在x0取極小值，所以

φ″(0)=(Ax,x)>0

所以A是正定矩陣。

反之，設(shè)A是正定矩陣，x0Ax=b的解，即Ax0-b=0，則由式(5.2)可得

φ(λ)=J(x0)+

(Ax,x)

=φ(0)+(Ax,x)>φ(0)

其中λ≠0，x≠0。

即J(x)在x0取極小值。由此可得到如下結(jié)論：當(dāng)矩陣A對(duì)稱正定時(shí)，二次函數(shù)J(x)的極小值和線性代數(shù)方程組Ax=b的解是一樣的。

J(x)是定義在空間Rn上的二次函數(shù)，稱為Rn上的二次泛函。5.1.2邊值問題轉(zhuǎn)化為變分問題

本章所討論的兩點(diǎn)邊值問題為下列形式：

(5.3)其中p(x)>0且是［a,b］上一次連續(xù)可導(dǎo)函數(shù)，q(x)≥0且是［a,

b］上的連續(xù)函數(shù)，f(x)是(a,b)上的平方可積函數(shù)。

對(duì)于邊值問題式(5.3)的求法，有兩種思路：一種基于極小位能原理，另一種基于虛功原理。

1.基于極小位能原理

構(gòu)造二次泛函

其中

令

a(u,v)=

得

J(u)=

a(u,u)-(f,u)

設(shè)H為所有平方可積且一階導(dǎo)數(shù)平方可積的函數(shù)中滿足左邊值條件u(a)=0的函數(shù)組成的集合，則邊值問題式(5.3)的變分問題為：求u*∈H，使J(u*)=

J(u)。a(u,v)是雙線性的，稱為雙線性泛函，即

a(c1u1+c2u2,v)=c1a(u1,v)+c2a(u2,v)

a(u,c1v1+c2v2)=c1a(u,v1)+c2a(u,v2)

而且a(u,v)是對(duì)稱的，即a(u,v)=a(v,u)，而a(u,v)的對(duì)稱性是由微分算子L決定的。可以證明，微分算子L是對(duì)稱正定的。和前面一樣，我們通過實(shí)變量λ的函數(shù)φ(λ)考察J(u)的極小值問題。任取v∈H，考慮實(shí)變量λ的函數(shù)。

下面引入極小位能原理：

定理5.1(極小位能原理)

設(shè)f是(a,b)上的連續(xù)函數(shù)，u*是二次連續(xù)可導(dǎo)函數(shù)，如果u*是邊值問題式(5.3)的解，則u*使J(u)達(dá)到極小值；反之，若H中的二次連續(xù)可導(dǎo)函數(shù)u*使J(u)達(dá)到極小值，則u*是邊值問題式(5.3)的解。

證明

a(u*,v)-(f,v)=如果u*是邊值問題式(5.3)的解，則Lu*-f=0，u*′(b)=0，因此對(duì)任意v∈H,

φ′(0)=a(u*,v)-(f,v)=0

所以

J(u*+λv)=J(u*)+

a(v,v)>J(u*)

即u*使J(u)達(dá)到極小值。反之，若u*使J(u)達(dá)到極小值，則對(duì)任意v∈H,

φ′(0)=a(u*,v)-(f,v)=

(Lu*-f)vdx+p(b)u*′(b)v(b)=0

若u*′(b)=0，則

(Lu*-f)vdx=0

根據(jù)變分法基本原理

Lu*-f=0

定理得證。以上我們對(duì)二次連續(xù)可微解u*(又稱古典解)建立了邊值問題和變分問題的等價(jià)性，但應(yīng)用中常涉及非二次連續(xù)可微的曲線，甚至是在某些點(diǎn)上不可微的分段連續(xù)可微函數(shù)，即只有一階導(dǎo)數(shù)。在實(shí)際應(yīng)用中此類J(u)也有意義。因此變分問題允許非光滑解，我們稱之為邊值問題的廣義解。簡(jiǎn)言之，古典解就是二次連續(xù)可微的廣義解。

2.基于虛功原理

微分方程兩端同時(shí)乘以v，并在［a,b］區(qū)間積分，得

(Lu-f)vdx=

其中

將其代入上式，得

即

a(u,v)-(f,v)=0

(5.4)由此可以看出，當(dāng)u∈H二次連續(xù)可導(dǎo)，且v∈H時(shí)，上述變分形式可化為

a(u,v)-(f,v)=

(Lu-f)vdx+p(b)u′(b)v(b)

若u是邊值問題式(5.3)的解，則u滿足上述變分形式；反之，若對(duì)任意v∈H，u∈H滿足上述變分形式，則u是邊值問題式(5.3)的解。由此我們得出以下結(jié)論定理：

定理5.2(虛功原理)

設(shè)u是二次連續(xù)可導(dǎo)函數(shù)，則u是邊值問題式(5.3)的解的充要條件是，u∈H且對(duì)任意v∈H，滿足變分方程a(u,v)-(f,v)=0。因?yàn)樯鲜龆ɡ淼耐茖?dǎo)過程中，方程算子不需對(duì)稱，所以上述定理可直接推廣為：

定理5.3

設(shè)u是二次連續(xù)可導(dǎo)函數(shù)，則u是邊值問題式(5.3)的解的充要條件是，u∈H，且對(duì)于任意v∈H，u滿足變分方程

a(u,v)-(f,v)=0

(5.5)

其中a(u,v)=。前面我們討論了如何將邊值問題轉(zhuǎn)化為變分問題，下面討論如何求解變分問題。一般情況下，無法求得變分問題的精確解，因此需要使用函數(shù)逼近，得出近似解。

我們知道任何函數(shù)都可精確表示為無窮多個(gè)基函數(shù)的線性組合的形式，在這里，我們用有限個(gè)基函數(shù)的線性組合近似精確解。設(shè)φ1，φ2，…，φn是一組基函數(shù)。由這些基函數(shù)的

線性組合可得到的函數(shù)集合記為Vn，其中某元素un表示為

5.2瑞茲-伽略金方法瑞茲法的目標(biāo)是選取系數(shù)ci，使二次泛函J(un)取極小值

J(un)=

a(un,un)-(f,un)=

cicja(φi,φj)-

cj(f,φj)

這是c1，c2，…，cn的二次函數(shù)，其中a(φi,φj)=a(φj,φi)。令

得線性代數(shù)方程組

a(φi,φj)ci=(f,φj),j=1,2,…,n

(5.6)求出ci，就得到近似解un。

伽略金方法是要求系數(shù)ci使un對(duì)任意v∈Vn滿足

a(un,v)-(f,v)=0或

a(φi,φj)ci=(f,φj),j=1,2,…

因?yàn)槿鹌澐ê唾ぢ越鸱▽?dǎo)出的線性方程組相同，因此稱該方程組為瑞茲-伽略金方法。需要指出的是，瑞茲法需要先把邊值問題化為等價(jià)的變分問題，而伽略金法不依賴于變分問題，所以當(dāng)且僅當(dāng)a(u,v)對(duì)稱正定時(shí)，瑞茲法才適用，對(duì)于a(u,v)非對(duì)稱正定的情況，只能采用伽略金法。

例5.1

用瑞茲-伽略金方法求解邊值問題

u″+u=-x,0<x<1

u(0)=0,u(1)=0

解取φi(x)=x(b-x)xi-1，將un表示成

un=

ciφi=x(1-x)(c1+c2x+…+cnxn-1)

令n=1，得u1=c1x(1-x)，由變分問題得

c1

(

+φ1)φ1dx=

-x2(1-x)dx可得

c1=

,u1=

x(1-x)

再令n=2，得u2=c1x(1-x)+c2x2(1-x)，由變分問題得

解之可得c1=，c2=，因此得到

該邊值問題的精確解為

表5-1列出了u1，u2，u*在x=，，的函數(shù)值。表5-1不同情況下u1、u2、u*的值注：根據(jù)右端常數(shù)項(xiàng)可以將邊值條件分為齊次邊值條件和非齊次邊值條件，如u(a)=0為齊次邊值條件，u(a)=d，d≠0為非齊次邊值條件；根據(jù)左端的形式可以將邊值條件分為第一類邊值條件、第二類邊值條件和第三類邊值條件，如u(a)=d為第一類邊值條件，u′(a)=d為第二類邊值條件,αu(a)+βu′(a)=d為第三類邊值條件，相應(yīng)的邊值問題稱為第一邊值問題、第二邊值問題和第三邊值問題。其中第一邊值條件又稱本質(zhì)邊值條件，第二邊值條件和第三邊值條件又稱自然邊值條件。對(duì)于非齊次本質(zhì)邊值條件u(a)=d1，可設(shè)

un(x)=d1φ0(x)+

ciφi(x)其中φ0(x)是滿足φ0(a)=1的任意函數(shù)，相應(yīng)的瑞茲-伽略金方程變?yōu)?/p>

cia(φi,φj)=(f,φj)-a(d1φ0,φj)

(5.7)

對(duì)于非齊次的第二邊值條件u′(b)=d2，則只需保留極小位能原理證明過程中的常數(shù)項(xiàng)p(b)u′(b)u(b)，相應(yīng)的瑞茲-伽略金方程變?yōu)?/p>

cia(φi,φj)=(f,φj)+p(b)d2φj(b)

(5.8)對(duì)于第三類邊值條件，為方便計(jì)算，設(shè)形式為

u(a)=d1

αu(b)+p(b)u′(b)=d2

以任意v(x)∈H乘以方程兩端，并對(duì)x在(a,b)上積分，利用分部積分得到

-p(x)u′(x)v(x)

p(x)u′(x)v′(x)dx+

q(x)u(x)v(x)dx=

f(x)v(x)dx

將右邊值條件兩端同乘以v(b)，得

p(b)u′(b)v(b)=d2v(b)-αu(b)v(b)代入上式，得

p(x)u′(x)v′(x)dx+

q(x)u(x)v(x)dx+αu(b)v(b)

=

f(x)v(x)dx+d2v(b)

即

a(u,v)=(f,v)+d2v(b)-αu(b)v(b)

相應(yīng)的瑞茲-伽略金方程為

cia(φi,φj)=(f,φj)+d2φj(b)-α

φi(b)φj(b)

(5.9)對(duì)于下列第三邊值條件所對(duì)應(yīng)的第三邊值問題，相應(yīng)的變分問題的二次泛函和瑞茲-伽略金方程的推導(dǎo)，作為課后習(xí)題，由學(xué)生們自行完成：

α1u(a)-β1u′(a)=d1

α2u(b)+β2u′(b)=d2

其中α1,α2≥0且不同時(shí)為零，β1,β2>0。上述瑞茲-伽略金方法將求解微分方程問題轉(zhuǎn)化成了求解線性方程組問題。隨著基函數(shù)的個(gè)數(shù)的增加，近似解un(x)也不斷逼近精確解u*(x)，當(dāng)n趨于無窮大時(shí)，所得到的解就是精確解u*(x)。利用瑞茲-伽略金方法得到的線性代數(shù)方程組的系數(shù)矩陣是一個(gè)稠密而且病態(tài)的矩陣，求解過程所需要的計(jì)算量和所占用的內(nèi)存空間相當(dāng)大，因此難以使用數(shù)值方法求解。5.3有限元法另外，以多項(xiàng)式基為例，由插值方法可知，當(dāng)基函數(shù)較多，多項(xiàng)式的次數(shù)較高時(shí)，龍格現(xiàn)象明顯，從而難以很好地逼近精確解u*(x)，于是產(chǎn)生很大誤差。為了克服這些缺點(diǎn)，

可以選擇支撐域足夠小的一系列函數(shù)作為基函數(shù)，而且只有相鄰兩個(gè)基函數(shù)的支撐域有重合，這樣由積分得到的線性方程組的大部分元素為零，也就可以解決上面的問題了。這就是我們下面將要介紹的有限單元法，簡(jiǎn)稱有限元法。本書僅討論求解一維問題的線性元。考慮兩點(diǎn)邊值問題

根據(jù)前面介紹的瑞茲-伽略金法，在變分過程中，可以從瑞茲法和伽略金法兩種方法得出上式的有限元方程。5.3.1從瑞茲法出發(fā)進(jìn)行推導(dǎo)

首先對(duì)求解區(qū)間進(jìn)行網(wǎng)格剖分，各節(jié)點(diǎn)為

a=x0<x1<…<xi<…<xn=b

每個(gè)小區(qū)間Ii=［xi-1,xi］稱為單元，單元長(zhǎng)度為hi=xi-xi-1

。最簡(jiǎn)單的情況，可用［a,b］上的一次分段函數(shù)近似微分方程的解，該一次分段函數(shù)由節(jié)點(diǎn)上的一組值u0=0,u1,u2,…,un通過線性插值公式

uh(x)=

ui-1+

ui,x∈Ii,i=1,2,…,n

得到,如圖5-1所示。圖5-1分段線性插值為了方便以后的計(jì)算，通常用仿射變換

ξ=Fi(x)=

將Ii變?yōu)棣屋S上的參考單元［0,1］。令

N0(ξ)=1-ξ,N1(ξ)=ξ

得

uh(x)=N0(ξ)ui-1+N1(ξ)ui,x∈Ii,i=1,2,…,n將上式代入泛函

J(u)=

(pu′2+qu2-2fu)dx

得

利用仿射變換，得

J(uh)=令

=aj-1,juj-1+ajjuj+aj+1,juj+1-bj=0

其中

(5.10)5.3.2從伽略金法出發(fā)進(jìn)行推導(dǎo)

在單元Ii，Ii+1對(duì)每一個(gè)節(jié)點(diǎn)構(gòu)造山形函數(shù)：

φi(x)=

φn(x)=

山形函數(shù)如圖5-2所示。圖5-2山形函數(shù)顯然φ1(x)，φ2(x)，…，φn(x)線性無關(guān)，可以用來作為逼近微分方程的解的基。uh可以表示為

uh(x)=

uiφi(x)

前面曾提到，與邊值問題相應(yīng)的雙線性泛函為

a(u,v)=

(pu′v′+quv)dx因此伽略金方程為

a(φi,φj)ui=

fφjdx,i=1,2,…,n

若借助仿射變換

N0(ξ)=1-ξ,N1(ξ)=ξ則對(duì)i=1,2,…,n-1

φi(x)=

顯然系數(shù)矩陣是三對(duì)角矩陣，即其中j=2,3,…,n-1。

右端項(xiàng)

上述三對(duì)角矩陣［aij］稱為總剛度矩陣，常數(shù)向bi稱為總載荷向量。

從上面的推導(dǎo)過程可以看出，伽略金法和瑞茲法得出了相同的有限元方程，但是通過伽略金法得出有限元方程更加方便，而且，由于伽略金法的適用范圍更廣，致使由其推導(dǎo)出

的有限元法的適用范圍也更加廣泛。

在前面介紹的有限元方法中，基函數(shù)由所有連續(xù)的折線函數(shù)組成，即由分段線性函數(shù)的線性組合構(gòu)成，該種方法稱為一次有限元，又稱線性元。要提高有限元方法的精度可以采用

更細(xì)密的區(qū)間分段(剖分)使單元長(zhǎng)度變小，節(jié)點(diǎn)參數(shù)增加；除此之外還有一個(gè)辦法就是采用更高次數(shù)的插值多項(xiàng)式作為基函數(shù)，這就是高次元。引進(jìn)了高次元的有限元法本書不再詳細(xì)論述，感興趣的讀者可以參考有限元方面的教材和專著。與其他方法相比，有限元法具有以下特點(diǎn)：

有限元法是以變分法為基礎(chǔ)的，是瑞茲-伽略金方法的進(jìn)一步發(fā)展，是一種很好的解決復(fù)雜區(qū)域和復(fù)雜邊值條件微分方程邊值問題的離散化方法；而且有限元法的各個(gè)計(jì)算環(huán)節(jié)都

便于標(biāo)準(zhǔn)化。對(duì)同一類型問題，不論幾何形狀和物理參數(shù)如何，以及采用不同的插值函數(shù)，都可用同一標(biāo)準(zhǔn)程序計(jì)算，這樣可以大大縮短解題周期。在工程計(jì)算中，并不是由前面所講的方式直接得到有限元方程，而是分析每個(gè)單元的局部二次型和單元矩陣，得到單元?jiǎng)偠染仃嚭洼d荷向量，再由單元?jiǎng)偠染仃嚭铣煽倓偠染仃?，?/p>

載荷向量合成總體載荷向量。這樣做的好處是易于用程序?qū)崿F(xiàn)。5.4實(shí)例介紹對(duì)比從上節(jié)得到的有限元方程可看出，ajj由

兩部分組成，bj由f(x)φj(x)dx和f(x)φj(x)dx兩部分組成，分別是相鄰兩個(gè)單元上的積分，而且有相同的形式，因此我們可以將總剛度矩陣和總載荷向量寫成單元?jiǎng)偠染仃嚭团c單元載荷向量和的形式。和單元?jiǎng)偠染仃嚨男问綖?/p>

其中單元載荷向量的形式為

其中可得有限元方程為

Ku=F

如果第一邊值條件非齊次，則需要對(duì)總剛度矩陣和總載荷向量做約束處理。如u(a)=d，可將單位載荷向量修改為

總載荷向量為

F=［d

…

］T

將總剛度矩陣修改為

K=

如果第二邊值條件非齊次，如

u′(b)=d

由相應(yīng)的瑞茲-伽略金方程可以看出，只需要將最后一個(gè)單位載荷向量修改為

=+p(b)d

即可。

對(duì)于非齊次的第三邊值條件

αu(b)+p(b)u′(b)=d

例5.2

用有限元法解下列二階常微分方程邊值問題

解顯然，在本例當(dāng)中

p(x)=

,q(x)=0,

f(x)=1

若將求解區(qū)間等分為四個(gè)單元，則節(jié)點(diǎn)坐標(biāo)為

x0=1,x1=

,x2=

,x3=

,x4=2每個(gè)單元長(zhǎng)度均為h=，下面兩式為山形函數(shù)系數(shù)

Ni(x)=

(xi+1-x)=4(xi+1-x)

Mi(x)=

(x-xi)=4(x-xi)

下面計(jì)算單元?jiǎng)偠染仃嚭洼d荷向量。單元?jiǎng)偠认禂?shù)為

單元?jiǎng)偠染仃嚍?/p>

單元載荷向量為總剛度矩陣為去掉第一行第一列，得總載荷向量為

F=

去掉第一個(gè)分量，得

F=得到有限元方程