華中科技大學(xué)計算方法課件-第七章常微分方程邊值問題數(shù)值解法

第七章常微分方程邊值問題數(shù)值解法一個由常微分方程及其邊界條件構(gòu)成的定解問題稱為常微分方程邊值問題，其常常出現(xiàn)在流體力學(xué)、材料力學(xué)、波動力學(xué)、電磁學(xué)及現(xiàn)代控制理論等學(xué)科中。一般而言，常微分方程邊值問題的解析求解是非常困難的。為了解決這類實際問題，人們通常采用數(shù)值解法。本節(jié)將主要介紹〃維常微分方程邊值問題fy\t)= ⑻,te(Q,b),tg(〃Q),貝6)=o ?)的打靶法、有限差分法及有限元法。

打靶法從形式上看邊值問題與初值問題之間的差別僅在于定解條件的不同，毫無疑問這兩類問題有著緊密的關(guān)聯(lián)。事實上，基于同一微分方程的邊值問題與初值問題有著共同的通解形式。因此，借助于初值問題數(shù)值算法來構(gòu)造邊值問題的數(shù)值算法是一個很直接的思路，而前者有很多有效的數(shù)值算法可供選擇。打靶法的基本思想是從不同的初值出發(fā)用已知的初值問題數(shù)值算法解相同的微分方程直到解的軌跡“打”到給定的右邊界點。為簡單起見，本節(jié)我們僅考慮如下7?維線性邊值問題/，=A（力加+。⑴，力e（電明\BqU?+Bby（b）=斗 （）其中4為和為為72X72的矩陣,q和〃為72維向量，而所用方法也可拓展到非線性邊值問題。

簡單打靶法對于一個常微分方程而言，不同的初值條件對應(yīng)于不同的解的軌線,即解是依賴于初始條件的。令此。)為方程/⑻=聞滿足初值條件沙(Q)=C的解,若貝"C)同時也滿足邊值問題(2)的邊界條件，那么v(1;c)正好就是邊值問題(2)的解。簡單打靶法要做的工作正是如何尋求這樣的初值。，使得Bac+Bby(b-c)=4 (3)根據(jù)常微分方程理論，我們知道線性邊值問題⑵的解可表示為y(t)=Y(t)c+v(t),t6(a,6), (4)其中y⑴，。⑴分別為基解矩陣及特解，cew1為一參數(shù)向量，且滿足TOC\o"1-5"\h\zY'=41)y,t€(q,b),Y(q)=/; (5)v'=A(t)v+q(t),v(a)=a, (6)這里/為〃級單位矩陣，acir通常可取為零向量。通過計算初值問題(5)-(6),我們可以得到邊值問題的解的表達式(4)。為了確定參數(shù)向量。，將(4)代入問題(2)的邊界條件可得Qc=力， (7)其中Q=Ba+場y⑹,3=1—Bav(a)-BbV(b).若邊值問題⑵存在唯一解，則Q是非奇異的，因此c是存在唯一的。當(dāng)。為零向量時，。正好為邊值問題的初值條件，在程序中算出c后，我們再次計算一次初值問題得到原邊值問題的數(shù)值解。

例7.1試分別取入=1,7,15,30,應(yīng)用簡單打靶法求解［0,1］上的邊值問題1oA21oA20(7r3+A27r)sin(7rZ)+(2A3+2Att2)cos(tt力)100\ /000\ /(e-2A+2e-A+3)/(e-A+2)000y⑼+1007/(1)= 00007 \010/ \(A(3-e-A))/(e-A+2)要求畫出數(shù)值解和精確解的對比圖，方程精確解為沙⑴=(〃1),〃，■),”"(t))T,其中u(t)=gA(t-l)_|_g2A(^-l)+g-At2+e-x+cos(M?

多重打靶法如果初值問題的積分區(qū)間較長，簡單打靶法的計算效果將會降I低。為克服簡單打靶法的這一缺點，下面我們介紹多重打靶法，其基本思想是首先將求解區(qū)間作網(wǎng)格分劃：q=力< <???<[N_i<In=A然后在每個小區(qū)間也一1,切上構(gòu)造微分方程的逼近解，最后由解的連續(xù)性及邊界條件可將這些數(shù)值解修補在一起得到一個全局解?？紤]線性問題(2),在區(qū)間也_1%]上，邊值問題的解可以表示為Ci.1)=匕⑴c1+a⑴，i=12…，N,(8)其中匕⑴，?、欧謩e為基解矩陣及特解，其滿足丫；=A⑴匕，匕(力―)=7, (9)耳=A(f)vi+q⑴，v他_口=0. (10)計算初值問題(9)-(10)得區(qū)間質(zhì)t,用上邊值問題的解的表達式⑻。為了確定參數(shù)向量c:=(%才,…將(8)代入問題(2)的邊界條件得TOC\o"1-5"\h\zCi—Y1gG_1= i<i<N-1, (11)Baco+BiXn⑸cn-i=〃一4Wv(6). (12)(11)-(12)可寫成如下線性系統(tǒng)Lc=r, (13)其中/-H(11)I \—K(12)IL- '1. ,?. ,-y?v-1(力n-1) i\Ba BbYN(b)Jr=(力(11)-2([2),N—i&n—i),〃一BbvN(b'))T.解線性系統(tǒng)(13)即得網(wǎng)格點力,上的數(shù)值解%=&。多重打靶法求解線性邊值問題的計算程序如下，其中輸入量表示意義與簡單打靶法基本相同，差別僅在力mes/z為求解區(qū)間上的全體網(wǎng)格點。算法7.2多重"靶醫(yī) [Y]=multishooting(basisfun,partifun,Ba,Bb,tmesh,eta)n=length(Ba);I=eye(n);T=zeros(n);N=length(tmesh)—1;L=zeros(N*n);r=zeros(N*n,1);Y=zeros(N+1,n);fork=1:N[t,y]=ode45(partifun,[tmesh(k),tmesh(k+l)],zeros(n,1));r(n*(k-l)+l:n*k)=y(Iength(t),:));forj=1:n[t,yJ=ode45(basisfun,[tmesh(k),tmesh(k+l)],I(：,j))；T(:,j)=y(length(t),:)5;endL(n*(k—l)+l:n*k,n*(k—1)+1:n*k)=-T;ifk<NL(n*(k—l)+l:n*k,n*k+l:n*(k+l))=I;endendL(n*(N-l)+l:n*N,1:n)=Ba;L(n*(N-1)+1:n*N,n*(N-1)+1:n*N)=Bb*T;v=r(n*(N-1)+1:n*N);r(n*(N—1)+1:n*N)=eta—Bb*v;c=L\r;fori=1:NY(i,:)=c(n*(i-1)+1:n*i)';endY(N+1,:)=(T*c(n*(N-l)+l:n*N)+v)\

-0.5將多重打靶法應(yīng)用于例7.1-0.5將多重打靶法應(yīng)用于例7.1中X=30的情形，可計算出相應(yīng)的數(shù)值解（見圖7.2）。由圖7.2可見，在計算精度方面，多重打靶法比簡單打靶法具有明顯的優(yōu)勢。t圖7.2例7.1（》=30）的精確解（虛線）及其多重打靶解（點劃線）.§7.2有限差分法有限差分法也是一種基于初值問題數(shù)值方法的邊值問題方法，但其不同于打靶法，有限差分法在利用初值問題數(shù)值方法離散邊值問題后，將求解區(qū)間上所有網(wǎng)格點處的數(shù)值解作為離散格式的未知量，從而通過求解表征該格式的線性或非線性代數(shù)方程組而獲得邊值問題的數(shù)值解。在本節(jié)，我們?nèi)匀豢紤]一階微分方程系統(tǒng)的邊值問題(1)-(2)。為了獲得邊值問題的系列離散格式，我們對求解區(qū)間[q㈤作如下剖分：n：Q=±0<±1<±2???<[N_l< "且記步長hj=ti—ti-i,%為V(切的逼近值?！?2]中點法與梯形法在親節(jié)，‘筑門先考高笥單的有限差分格式，即中點法和梯形法。分別應(yīng)用中點法和梯形法到非線性邊值問題(1)-(2),我們依次可得其問題的離散格式=f,1/2,；+4一1))('=1,2, g?o,VN)=0" (14)仇=才/9'納)+/9—1，%—1)]&=1,2,...,N),g(go,Vn)=0." (15)上述計算格式可視為以向量)=(詔，端，…，川)丁為未知量的“(N+1)維非線性方程組，其可用Newton迭代法求解。為敘述簡單見，這里我們僅給出梯形法的Newton迭代格式，中點法的Newton迭代格式類似可得。記NhU：="J—~ 仇)+/(力:—i,%—1)],i=l,2,…，NF⑼二(M-,9必，…,明曲g(yo,—則中點法的Newton迭代格式為尸?〃'）￡=-W獷，+1=/〃+您加=o,1,2,???,W（16）其中尸（。為9（?）的Jacobi矩陣，m為迭代次數(shù)。若引入記號：￡=以…，雜,4地=曲，，，），氏=&（琢,-）,Bb=必端淵.則迭代格式（16）可進一步寫成f 幻&+4力-1）￡-1]=—Nw八2=1,2,…，N,]&a'+場就=—。（婚,歿）；、次+1=次+。￡=0,…，N.（17）在格式（⑺的每步迭代中，通過前二組方程求出&然后將其代入第三組方程求出當(dāng)前逼近值必+1 =[（鱷+丫,破+1）丁,…,（第+1）平。其計算程序如下：算法7.3中點方法[y]=midpoint(rfun,prfun,gfun,pgfun1,pgfun2,tmesh,yO)[N,n]=size(yO)N=N-1;I=eye(n);L=zeros(n*(N+l));r=zeros(n*(N+l),l);tol=1.d—12;xi=ones(n*(N+l),1);y=yO;while(norm(xi)>tol)forj=1:Nhj=tmesh(j+l)—tmesh(j);A=feval(prfun,tmesh(j)+hj/2,(y(j,:)'+y(j+l,:)*)/2);L(n*(j-1)+1:n*j,n*(j—1)+1:n*j)=-1/hj*1—A/2;L(n*(j-l)+l:n*i,n*j+l:n*(j+1))=1/hi*I-A/2;r(n*(j-l)+l:n*j)=-(y(j+1,：)'-y(j,:)')/hj+...feval(rfun,tmesh(j)+hj/2,(y(j,:)'+y(J+1,:)')/2);endL(n+N+l:n*(N+l),l:n)=feval(pgfunl,y(l,:)*,y(N+l,:)));L(n*N+l:n*(N+l),n*N+l:n*(N+1))=feval(pgfun2,y(l,y(N+1,:)))；r(n*N+1:n*(N+l))=—feval(gfun,y(1,:)',y(N+l,:)');xi=L\r;fork=1:N+1y(k,:)=y(k,:)+xi(n*(k—l)+l:n*k)';endendend上述程序中，輸入量Tfun和prfun分別表示函數(shù)三祟;gfun,pgfuni,pgfun2分別表示函數(shù)縹署,寫21mes/z表示求解區(qū)間上所有網(wǎng)格點［如由，…，6］，y。表示網(wǎng)格點上的初始迭代值。

取步長九=1/(10*%)(，=0,1,2,3),初始迭代值為零，分別應(yīng)用中點法和梯形法求解該邊值問題可計算出其數(shù)值解的全局誤差及收斂階(見表7.1),其中全局誤差及收斂階的表達式依次為、 ?/、 , 1 「err(h)'err(h):=maxkzL—uA,p:=logo —.-o<i<N1kl)上? 82Lerr(/z/2)Jh中點法梯形法err(h)Perr(h)P1/106.3341e-4-5.8442e-4-1/201.5938e-41.99071.4632e-41.99791/403.9867e-51.99923.6594e-51.99941/809.9698e-61.99969.1516e-61.9995表7.1中點法和梯形法解例7.2的數(shù)值結(jié)果.由表7.1可見，中點法和梯形法均只有2階精度，并且二者的計算效果相近。

§7.2.2Runge-Kutta方法由上可知，中點法和梯形法是二種低精度方法。為獲得高精度的邊值問題方法，我們考慮修改求解初值問題的Runge-Kutta方法，使其適用于邊值問題(1)-(2)。將Runge-Kutta方法應(yīng)用于邊值問題⑴可得NfJij:=fij-玲,m—i-hif=0,<i<N,1<j<s,NhUi：=Vi-Ui—i-hitbjfij=0,1<i<Ng(y^yN)=0,、 j=i(18)這里%=tji+hlCjo該方程組中的未知向量y：=(yLfLy:4,…，氐武)「力:=(危忌…，息),可通過Newton迭代法獲得，其計算格式為F\ym^=-F(r),ym+1=r+€,m=0,l,2,... (19)

其中尸(?)為其中尸(?)為9(。的Jacobi矩陣，皿為迭代次數(shù)，且r,， r,， r,，-,NhfN,NhyN,g(y^yN))&&+場Sn=—。（鱷,第），次+1=次+&,0<z<7V,f/1= +%l<i<N.則迭代格式(19)可寫成￡=(8,端則迭代格式(19)可寫成既然(20)中的第一個方程可寫為口=W”￡—1+叱TQz,USN,因此將其代入(20)中的第二個方程可得金一1總—1=/八1</<A^,其中八 『2=/+h,DW^Vi,n=hiDW^Qi.從而迭代格式(20)可寫成'a一aWi-i=/逅1<z<a7-,,&&+36v=-。(鱷,第)， 。、\少+1=次+&,0<z<7V, ⑷)、犁+1=花+叱-1匕￡_]+嗎-依，l<i<N.在格式(21)的每次迭代中，先通過前兩個方程計算出8然后代入后兩個方程完成一次迭代，其計算程序見算法7.4。在該算法中，輸入量基本與算法7.3的同名輸入量相同，唯一不同的是這里多了一個輸入量加，它表示右端函數(shù)在整個區(qū)間上級點處的初始迭代值(足…，后了。

§7.2.3線性多步法本節(jié)，我們考慮修改求解初值問題的線性多步法，使其可應(yīng)用于邊值問題(1)-(2)。邊值問題線性多步法的基本思想是：首先在求解區(qū)間二端及其附近分別給出島個和松個逼近值，然后在求解區(qū)間中部的每個網(wǎng)格點右處用線性多步法離散，將關(guān)于逼近解｛功｝的方程組與邊界方程=0聯(lián)立，通過適當(dāng)方法求解該聯(lián)立方程組而獲得原邊值問題(1)-(2)的數(shù)值解。邊值問題線性多步法的一般形式為):Qj+kiVi+j=h〉:Bj+kJ&i+j,yi+j)，￡=瓦,?..，N一例，(22)j=~ki j=~ki其中ki+k2=k,to=a,ti= +h,h=(b-a)/N,僅為g(幻的逼近。為了求解方程組(22),我們?nèi)孕枰猭個逼近值比,陰，…，班1-1，以-比+1，…，以，其可由k-1個與(22)同階的線性多步公式TOC\o"1-5"\h\zk—i k—i￡ = 力比+j),2=1,... -1, (23)j=~i j=~ik—i k—i):a；［jVN-k+i+j=h): jf(tN—k+i+j'>UN—k+i+j),〃=N—%+l,???,Nj=~i j=~i(24)及邊界條件g("o,UN)=。 (25)給出。至此，我們即得到求解邊值問題(1)-(2)的線性多步法(22)-(25)。

為實現(xiàn)格式(22)-(25)的計算，我們將(22)-(24)寫成矩陣形式(A0l)y=h{B^I)F(y1 (26)其中V=(湍,對,必)。F⑼=(/(%,yo￥\f(ti,城,…Vn)t>t,a。(s-A:2+1)a。(s-A:2+1)“0a1a(s—極+i)(s-fc2+l)ak(s)%,。心）分別代之JJ4有相同形式，僅需將月中的以%徘）即可獲得。聯(lián)立（26）與%,。心）分別代之JJ4有相同形式，僅需將月中的下面，我們介紹幾個常用的邊值問題線性多步法，其推導(dǎo)及相關(guān)理論可參見Brugnano和Trigiante的專著。三階廣義向后微分公式(GBDF)(能為+i+3yL6%—1+yn-2)=6力，，=2,3,…，N—1,\5(一明+6g2-3陰-2t/o)=九九I—18^-1+§Vn-2—2gyv_3)=hf^.五階廣義Adams方法(GAM)i%一%-1=含(-19力―2+346%_1+456力―7"+i+1"+2),，=2,.「N-2加一班=含(一19人+106/3-264/2+646力+251/0),W-i-Un-2=盒(―19冊+346/ar—i+456/y-2—74/n―3+ll/v-4),、Un—yN-i=含(251/n+646/tv-i—264/n—2+106加-3—19/2v_4).

?四階擴展的第二類梯形方法(ETR2)(吉(%+1+9仇一9%_l仇—2)=亨(力+力一1),i=2,…)N—1,S貴(一為+9n2+9陰-17"))=1(3/i+用，［貴(VN-3—9yN-2—9yN-i+17〃n)=亨(3/n-i+/at).?六階方法(TOM)'看(11如1+27%一27比—1—llyi-2)=，(力+1+9fi+9/i-i+力一2),2=2,...,7V—1,<叭一處=儡(27右—173A+482/3-798/2+1427力+475/0),Vn-yN-i—i^(27/n-5—173加-4+4826-3—798/n-2、 +1427^-1+475/^).線性多步法求解非線性邊值問題(1)-(2)的計算程序見算法7.5,其中前七個輸入量的意義與算法7.3相同。算法7.5線性多步法function[y]=bvm(rfun,prfun,gfun,pgfun1,pgfun2,tmesh,yO,A,B)[N,n]=size(yO);N=N—1;I=eye(n);h=tmesh(2)—tmesh(1);L=zeros(n*(N+1));J=L;r=zeros(n*(N+l),l);tol=1.d-16;xi=ones(n*(N+l),1);y=yO;while(norm(xi)>tol)forj=1:N+1Y(n*(j-l)+l:n*j,l)=y(j;F(n*(j—l)+l:n*j,:)=feval(rfun,tmesh(j),y(jJ(n*(j-l)+l:n*j,n*(j-l)+l:n*j)=feval(prfun,tmesh(j),y(j,:)')endL(1:N*n,:)=kron(A,I)—h*kron(B,I)*J;r(1:N*n)=—kron(A,I)*Y+h*kron(B,I)*F;L(n*N+l:n*(N+l),l:n)=feval(pgfunl,y(l,:)5,y(N+l,:)');L(n*N+1:n*(N+1),n*N+1:n*(N+1))=feval(pgfun2,y(l',y(N+lr(n*N+1:n*(N+l))=—feval(gfun,y(l ',y(N+l,:)');xi=L\r;fork=1:N+1y(k,:)=y(k,:)+xi(n*(k-1)+1:n*k));endendend

為測試其算法的有效性，我們分別取步長h=1/(5*20(i=0,1,2),依次應(yīng)用上述四種計算格式于例7.2中的邊值問題，其數(shù)值解的全局誤差err(h):=maxerr(h):=maxOWigNI“閻~Ui\和收斂階irerr(h)一

P=10g2岸向列于表7.3中。由該計算結(jié)果可見，其算法關(guān)于非線性邊值問題是非常有效的。hGBDF-3ETR2-4GAM-5TOM-6err[h)perriji)perr(Ji)perr(Ji)p1/56.3674e-4-3.1464e-5-1.3510e-5-7.5878e-6-1/109.7062e-52.71372.0233e-63.95894.9753e-74.76314.4113e-87.42631/201.4396e-52.75321.2918e-73.96931.7827e-84.80269.7162e-105.5047表7.3線性多步法解例7.2的數(shù)值結(jié)果.§7.2Ritz-Galerkin方法源于力學(xué)中的變分原理有著重要的理論和實踐價值，變分原理主要有極小位能原理和虛功原理。在特定的函數(shù)空間設(shè)置下，微分方程邊值問題可轉(zhuǎn)化為與之等價的變分問題，通過求解相應(yīng)的變分問題，我們可獲得原微分方程邊值問題的解。其中，基于極小位能原理產(chǎn)生的方法稱為Ritz方法，而基于虛功原理產(chǎn)生的方法則稱為Galerkin方法，有時人們也籠統(tǒng)地稱這二種方法為Ritz-Galerkin方法。抽象空間為了方便本節(jié)及后面第九章的應(yīng)用，這里我們首先介紹一些抽象空間概念和常用的積分不等式。設(shè)有區(qū)域。ursa。為其光滑邊界，C=為其閉包，則我們可給出如下抽象空間概念及其常用結(jié)論：C『(0)空間：在。上無窮次可微，在3。的某鄰域上恒為零的全體函數(shù)。C，(。)空間：在Q上m次連續(xù)可微，在3。的某鄰域上恒為零的全體函數(shù)。V(Q)空間(1<p<oo)：在。上p次絕對可積的全體函數(shù)。易知，空間117(。)={/:11/11^(0)<oo}是一個Banach空間。當(dāng)p=2時，L2(。)空間中的內(nèi)積和范數(shù)可分別定義為（工。）=[fgdo,||/||=\[\f\2d（J

Jo VJo此外，在LP(O)空間中，下列不等式成立：Minkowski木琴式：設(shè)/,gGLP(Q),則11/+^||v(Q)<||/||LP(n)+II^IIlp(Q)-Holder不等式：設(shè)/e叫0),ge*0),+1=1,貝lj%eL!(Q),且llfgllu⑼，II/IIlp(q)I|5,IIl9(q)-當(dāng)P=q=2時，Holder不等式即為Cauchy-Schwarz不等式II/^IIv(q)<||/IIl2(q)||sf||l2(q).變分法基本引理：設(shè)/eL2(Q),且有If(x)^[x}d(y=0,V@eC^°(Q),Jo則/(力)在。的閉包。上幾乎處處為零；特別，若/eC(C),且使得上式成立，則/Q)三0(16。)。一階廣義導(dǎo)數(shù)：設(shè)/eL2(Q),若存在g6L2(Q),使得[g(a)@⑶dcr=-/ V3GC^°(Q),Jo JoOXi其中i …，2n)丁，則稱g是/關(guān)于◎的一階廣義導(dǎo)數(shù)，且記為小/儂)。同階廣義導(dǎo)數(shù):設(shè)/GL2(Q),若存在qGL2(Q),使得Ig(N)9Q)dcr=(-1)，jf(x)da^d(r,V砂.玲⑼,Jo Jo其中7=(11,11,???,7冗)了，且q=(Q1,02,…，an)€RM,|a|=￡%,3°融=7r為一百FuX-\???uXjiJ=1 1則稱g為了的同階廣義導(dǎo)數(shù)，且記為V/(c)。顯然廣義導(dǎo)數(shù)有別于經(jīng)典導(dǎo)數(shù)，前者描述的是全局性質(zhì)，而后者描述的是局部性質(zhì)。二者之間有如下關(guān)系：若/eC悶(Q),則f的廣義導(dǎo)數(shù)與經(jīng)典導(dǎo)數(shù)相同。

Ritz方法考慮一維橢圓邊值問題fd/Ly：=一丁加(1)“(叫+q(x)y(x)=f⑺,dxieC=(Q,b), (27)y⑷=o,y'(b)=o,其中pq6CX(Q),p(x)>p():=minp(x)>0,q(x)>0(reGQ),f6L2(Q)o對于滿足上述邊界條件且在。上二次連續(xù)可微的任意兩函數(shù)次I),或乃，利用分部積分公式可得(Ly,z)=/[―而(4)+磯2五=/(py'z'+qyz)dx.

Ja Ja記上式右端為Q(/N),即有(Ly、z)=a(y,z). (28)泛函q?,z)將在我們建立邊值問題(27)的變分原理時起著舉足輕重的作用，特別泛函Q(〃N)還具有如下重要性質(zhì)：1）線性性：若a,0為任意實數(shù)，則Q（W/1+By*z）=QQ（91,z）+戶。僅2,z）,Vgi,1/2,Ze回1（。）,

Q（V,azi+戶Z2）=aa（/2i）+Pq（/22）,%,孫3G皿］。）2）對稱性：a（y,z）=a（2,y）,\fy,zWH'0）3）連續(xù)性：存在一個常數(shù)A/>0,使得|q?,z）|<mi訓(xùn)111Hli,Vy,zg明。）4）強制性：存在常數(shù)7>0使得Q（-y）>y\\y\\tVge唯（。）這里={y:yG?。ā＃恆）=0}。顯然Hj（Q）CIHI^（Q）crf（Q）.

強制性的證明：口/eHX。）有因此,由Cauchy-Schwarz不等式得1fb ]正F，0⑵|2七+，b\y\x)\2dx>刊Ml；,這里，=minG，島7｝。故證明VzeHlX。），我們有a(n*,N)-(/,z)= (py'*z'+qy*z-fz)dxJa,b fb d , -i=(py'*z)+ -丁(01)+qy*—于zdx(30)QJaLdx 」=P(b)&(b)z(b)+[(Ly*-f)zdx\Ja若科是邊值問題(27)的解，則有Ly*-于=0,y'*(b)=0.從而，據(jù)(30)得a(g*,z)-(￡z)=0,\/ze (31)由上式及q⑵z)的強制性可知：VA€肽\{0}及VzGHX0\{O},泛函J?)滿足J(y*+1?)=J(U*)+z)>J(y*). (32)因此，J(U*)=ruinJ(y)。y典⑼

反之，若J(g*)=minJ(y),則當(dāng)記①(八)：=J(g*+M)= +A[a(?/,,z)-(/,z)]+—a(z,z),AeR時，由(30)我們有rbTOC\o"1-5"\h\z①/(0)=p(b)山⑼2(b)+ (Ly*—f)zdx=0,Vzs (33)且因此 ?/{Ly,-f)zdx=0,\/zGCg0(Q). (34)Ja應(yīng)用變分法基本引理于(34)得Ly*-f三0,VxeQ. (35)將(35)代回(33)得P(b)“*(b)z(b)=0,\/zeH^(Q).而p(b)>po>0>且取z(2)=x—q時有z伍)>Oo因此暝(b)=0o綜上所述，g*為邊值問題(27)的解?！鲈摱ɡ硪蠼馔舛芜B續(xù)可微，此時我們稱之為邊值問題(27)的經(jīng)典解。而就變分問題(29)而言，其解外不必滿足如此強的光滑條件，因此我們稱變分問題(29)的解為邊值問題(27)的廣義解或弱解。以下，我們基于極小位能原理，通過求解變分問題(29)而獲得原邊值問題(27)的逼近解。設(shè)Um是回X。)的加維子空間(稱之力試探函數(shù)空回)：｛⑵⑵｝陶是?！钡囊唤M基函數(shù)，則心中任一兀素％論(1)可表為m=丹 (36)i=l

今選取系數(shù)｛?！?，使J(加0取極小?？紤]到] m mJ(Vm)=5砌的Vm)~(7,Vm)=j 仍)。嗎一￡/(￡丹)，J=1 j=l是｛4｝的二次函數(shù)，且a(g,仍)具對稱性，則令dcjdcjo,J-1,2,???,m可得m￡qS,=(/,0)，j=12…9. (37)2=1解此關(guān)工｛&｝的線性方程組，且將其解代回(36)即得變分問題(29)的逼近解加，這就是Ritz方法。

§7.3.3Galerkin方^去由定理7.1的證明過程，我們可知：若ye?2(0)是邊值問題(27)的解,則有z)-(/,z)=0,\/zG (38)反之，若?2(。)是變分方程(38)的解，貝IJ由(30)得p[b}y'(b)z(b)+(Ly—f)zdx=0,\/zEJa利用該方程及定理7」后半部分的證明可推得v是邊值問題(27)的解。因此，我們有以下變分原理。定理7.2設(shè)函數(shù)vGC2(Q),貝切是邊值問題(27)的解當(dāng)且僅當(dāng)V是變分方程(38)的解。

由于變分方程(38)的左端在物理學(xué)中常常表示虛功，因此我們稱定理7.2為虛功原理，而其方程的解也稱為邊值問題(27)的廣義解或弱解。以下，我們基于虛功原理，通過求解變分方程(38)而獲得原邊值問題(27)的逼近解。設(shè)%是皿乂。)的“維子空間，｛處⑺｝黑1是的一組基函數(shù)，則。加中任一元素心⑶可表示為(36)。今選取系數(shù)｛踴｝,使得%〃(如滿足變分方程TOC\o"1-5"\h\zz)=(/,z),以eUm. (39)m記N=￡dj燈(Vd,GR),將該式及(36)代入(39)得，=im m m￡CiQ?,R)=￡4Cj=l i=l j=l由功的任意性得m￡一(-/)■=(力㈤,J=1,2,???,m.2=1該式即為(37),解該方程即可獲得(38)的逼近解而，此即為Galerkin方法。

Ritz方法與Galerkin方法從不同途徑獲得了方程(37),因此我們稱方程(37)為Ritz-Galerkin方程,而將Ritz方法與Galerkin方法統(tǒng)森為Ritz-Galerkin方法。但是，我們必須注意，Galerkin方法比Ritz方法具有更廣的可應(yīng)用范圍。如：考慮比(27)更一般的邊值問題(dLy：=一二加OWO)]+ry\x)+q(x)y(x)=f(x),dX -Q,b), (40)v(q)=o,y'(b)=o,其中r,p,q￡CX(Q),p(x)>^>o>0,q(x)>0(rreQ),/6L2(Q)O若記rba(y,z)=(py'z’+ry'z+qyz)dxyJa則我們可建立(40)的虛功原理，并由此導(dǎo)出其Galerkin方法。但是，由于此時的泛函a?,z)不再具有對稱性及強制性，因此不可能建立(40)的極小位能原理，故而也不可能導(dǎo)出其Ritz方法。

例7.4取基函數(shù)2(c)=sin(27ra;)(i=1,2,???,m),用Galerkin方法求解非齊次邊值問題(—y"(x)+y⑺=x2,xe(0,1),［貝0)=0,貝1)=L解令隊=%u=g—如，則欲解的非齊次邊值問題可化為齊次邊值問題(—u"(x)+u(x)=x2—x,xG(0,1),

(“(0)=0,M⑴=0.m設(shè)該問題的6次Galerkin逼近解為〃6(幻=￡ ⑶，則其滿i=l足Ritz-Galerkin方程m￡Q(一,—)G=(X2-x,(pj),j= ,mi=l

這里這里.2-x,(Pj)=J(x2—x)(pj(x)dx=晨)3[(-I)，-1].因此 _s,、大一c_(/_7,Q)_f-12+1]，，為奇數(shù)，' 。(內(nèi)白) [o,，為偶數(shù)、口 Im+1I—8sin[(2A;—l)7ra;]S=A[(2k—1閉3{[(27—1閉2+1}，這里用表示取整函數(shù)。故原邊值問題的Galerkin逼近解為—8sin[(2A;-1)ttx]

[(2k—1)汗血(2心一1用2+1}.§7.3.4有限元法邊值問題的Ritz-Galerkin方法存在著若干弊病，首先，當(dāng)形成Ritz-Galerkin方程時，我們需要計算大量的積分，例取加維試探函數(shù)空間/，則需要計算*加（加+3）個（對稱情形）到加（小+1）個俳對稱情形）積分；其次，若Ritz-Galerkin方程是通過取傳統(tǒng)基函數(shù)形成的，則該方程組往往是非稀疏的大型病態(tài)方程組，其求解十分困難。為克服其困難性，本節(jié)介紹有限元方法，該方法與Ritz-Galerkin方法主要區(qū)別在于其選取局部基函數(shù)的技巧。目前，有限元方法已被廣泛應(yīng)用于結(jié)構(gòu)力學(xué)、固體力學(xué)和流體力學(xué)等工程領(lǐng)域，其求解的大致步驟如下：Step1.將邊值問題轉(zhuǎn)化為變分問題;Step2.將求解域剖分成適當(dāng)形狀的單元;St叩3.構(gòu)造單元形狀函數(shù),形成有限元空間;Step4.形成有限元方程,并