第二章有限差分法主講人：胡才博中國科學院大學地球科學學院中國科學院計算地球動力學重點實驗室第二章有限差分法主講人：胡才博第二章有限差分法2.1有限差分法基礎2.2網(wǎng)格剖分2.3差分格式2.4差分方程2.5應用實例第二章有限差分法2.1有限差分法基礎1.地球內(nèi)部介質(zhì)，不僅存在縱向非均勻結構（一維地球模型），也存在橫向非均勻結構（不同塊體、斷層系統(tǒng)）；2.幾何模型也呈現(xiàn)出相當?shù)膹碗s性；3.另外，邊界條件和初始條件對于不同問題具有特殊性。解析方法的局限性Hu,C.,Y.Cai,andZ.Wang(2012),Effectsoflargehistoricalearthquakes,viscousrelaxation,andtectonicloadingonthe2008Wenchuanearthquake,JournalofGeophysicalResearch,117,B06410,doi:10.1029/2011JB009046.(SCI,IF:3.303)汶川大地震的動力學成因1.地球內(nèi)部介質(zhì)，不僅存在縱向非均勻結構（一維地球模型），對于存在復雜介質(zhì)和幾何、特殊邊界條件和初始條件的實際地質(zhì)問題，一般不存在解析解，需要近似的數(shù)值求解方法。有限差分方法是地球物理方法中最常見的一種。有限差分方法（FiniteDifferenceMethod,FDM）是計算機數(shù)值模擬最早采用的方法，至今仍被廣泛使用。有限差分方法的基本特點該方法是一種直接將微分問題變?yōu)榇鷶?shù)問題的近似數(shù)值解法，數(shù)學概念直觀，表達簡單，是發(fā)展較早且比較成熟的數(shù)值方法。2.1有限差分法基礎對于存在復雜介質(zhì)和幾何、特殊邊界條件和初始條件的實際地質(zhì)問題有限差分法以變量離散取值后對應的函數(shù)值來近似微分方程中獨立變量的連續(xù)取值。我們放棄了微分方程中獨立變量可以取連續(xù)值的特征，而關注獨立變量離散取值后對應的函數(shù)值。有限差分法的具體操作分為兩個部分：(1)用差分代替微分方程中的微分，將連續(xù)變化的變量離散化，從而得到差分方程組的數(shù)學形式；(2)求解差分方程組。有限差分方法的基本原理該方法將求解域劃分為差分網(wǎng)格，用有限個網(wǎng)格節(jié)點代替連續(xù)的求解域。有限差分方法以Taylor級數(shù)展開等方法，把控制方程中的導數(shù)用網(wǎng)格節(jié)點上的函數(shù)值的差商代替進行離散，從而建立以網(wǎng)格節(jié)點上的函數(shù)值為未知數(shù)的代數(shù)方程組。2.1有限差分法基礎有限差分方法的基本原理該方法將求解域劃分為差分網(wǎng)格，用有限個有限差分法的主要內(nèi)容建立地球物理問題的離散有限差分模型（1）如何根據(jù)問題的特點將定解區(qū)域做網(wǎng)格劃分；（2）如何在所有網(wǎng)格節(jié)點上用有限差分格式對導數(shù)求近似，對函數(shù)、初始條件和邊界條件求近似；（3）如何把原方程離散化為代數(shù)方程組，即有限差分方程組。2.從理論上研究有限差分模型的形態(tài)，以保證計算過程的可行性和計算結果的正確性（1）解的相容性；（2）解的穩(wěn)定性；（3）解的收斂性。3.如何數(shù)值求解差分方程組有限差分法的主要內(nèi)容建立地球物理問題的離散有限差分模型2.從網(wǎng)格剖分就是研究區(qū)域和邊界的離散化1.矩形分割2.三角形分割3.極網(wǎng)格分割2.2網(wǎng)格剖分網(wǎng)格剖分就是研究區(qū)域和邊界的離散化2.2網(wǎng)格剖分對地球物理問題的連續(xù)求解區(qū)域通過網(wǎng)格劃分離散為空間上得一系列網(wǎng)格點，接下來需要利用一定的差分格式對偏微分方程組中的導數(shù)用差商進行近似，從而將偏微分方程組離散化為差分方程組。對于函數(shù)f(x)，通常意義下的導數(shù)（微商）定義為：當dx→0時，以上三種形式都是微商的正確定義。如果dx是有限的，如何給出微商的近似定義？對地球物理問題的連續(xù)求解區(qū)域通過網(wǎng)格劃分離散為空間上得一系列2.3差分格式用Taylor級數(shù)展開可以給出微商的近似形式。對于連續(xù)函數(shù)f(x)，它在相鄰點上的值f(x+Δx)和f(x-Δx)可以用Taylor級數(shù)展開為dx變?yōu)棣如果Δx很小，f(x)可微，則以上級數(shù)收斂。次數(shù)越高，收斂級數(shù)的項的絕對值越小。由（1）得到，（1）（2）（3）（4）2.3差分格式用Taylor級數(shù)展開可以給出微商的近似形式式中的O(x)項表示忽略掉的所有項中的最大項的量級是Δx，也就是說，忽略掉這些項帶來的誤差中的最大項和Δx成正比。由（4）給出導數(shù)的一階精度（firstorderaccurate）近似為：（4）（5）（5）式稱為向前差分格式（forward-differenceformula）由（2）式得到由（7）式得到導數(shù)的另一個一階精度近似：（6）（7）（8）（8）式稱為向后差分形式（backward-differenceformula）。式中的O(x)項表示忽略掉的所有項中的最大項的量級是Δx，也（1）式減去（2）式，得到：（9）式中的O(Δx2)項表示忽略掉這些項帶來的誤差中的最大項和Δx2成正比。由（9）式得到導數(shù)的二階精度（secondorderaccurate）近似為：（10）式稱為中心差分形式（central-differenceformula）。（9）（10）（1）（2）（1）式減去（2）式，得到：（9）式中的O(Δx2)項表示忽（1）（2）（1）式和（2）式相加，得到：（12）式稱為二階導數(shù)的二階精度中心差分形式。忽略Δx的四次方及更高階項(11)(12)f(xi+h)-f(xi):節(jié)點xi的一階向前差分f(xi)-f(xi-h):節(jié)點xi的一階向后差分f(xi+h)-f(xi-h):節(jié)點xi的一階中心差分前后是相對x軸正方向而言（1）（2）（1）式和（2）式相加，得到：（12）式稱為二階總結：1、向前差分形式：2、向后差分形式：3、中心差分形式：單側，一階精度單側，一階精度對稱，二階精度對于二階導數(shù)二階精度對一階導數(shù)總結：1、向前差分形式：2、向后差分形式：3、中心差分形式：定解問題的有限差分解法1．離散x=ih,y=jh,i=0,±1,±2,….±n,h:步長（正方形的邊長）2．根據(jù)泰勒級數(shù)建立差商格式：對于一維情況:在x處的一階導數(shù)可以用3.建立和求解差商方程組3.建立和求解差商方程組差分格式的另一種推導為了尋求更精確的差分格式，我們引入兩個待定常數(shù)，由泰勒展開，構造如下關系式(i,j)(i+1,j)(i-1,j)01234(i,j)(i+1,j-1)(i-1,j-1)(i,j+1)(i+1,j+1)(i-1,j+1)i-1ii+1j-1jj+1h1h3h2h4差分格式的另一種推導為了尋求更精確的差分格式，我們引入兩個待為了尋求更精確的差分格式，我們引入兩個待定常數(shù)，由泰勒展開，構造如下關系式回代（1）中，舍去高階項一階偏導數(shù)中心差分的推導（1）為了尋求更精確的差分格式，我們引入兩個待定常數(shù)，由泰勒展開，（1）二階偏導數(shù)差分的推導回代（1）中，舍去高階項二階偏導數(shù)差分公式（1）二階偏導數(shù)差分的推導回代（1）中，舍去高階項二階偏導數(shù)一個例子：等步長一個例子：等步長一個例子：局部節(jié)點離散化方程總體節(jié)點離散化方程總體節(jié)點離散化方程f=0時，變?yōu)椴此煞匠蘤=q=0時，變?yōu)槔绽狗匠桃粋€例子：局部節(jié)點離散化方程總體節(jié)點離散化方程總體節(jié)點離散化（1）第一類邊界條件（a）直接轉移法在圖中網(wǎng)格是按正方形分割，步長為h。0點為靠近邊界G的一個網(wǎng)格節(jié)點，1和2為邊界節(jié)點。我們?nèi)∽羁拷?點的邊界節(jié)點1上的函數(shù)值作為0點的函數(shù)值。即取φ0≈φ1。這種方法稱為直接轉移法，又稱為零次插值法。

邊界條件的離散化的處理狄利克萊問題（1）第一類邊界條件（a）直接轉移法邊界條件的離散化的處（b）線性插值法先判斷x方向的邊界節(jié)點1和y方向的邊界節(jié)點2哪一個更靠近0點。如果1更靠近0點，則可以用x方向的線性插值給出0點的函數(shù)值如果2更靠近0點，則可以用x方向的線性插值給出0點的函數(shù)值（b）線性插值法如果1更靠近0點，則可以用x方向的線性插值（c）雙向插值法(i,j)(i+1,j)(i-1,j)01234(i,j)(i+1,j-1)(i-1,j-1)(i,j+1)(i+1,j+1)(i-1,j+1)i-1ii+1j-1jj+1h1h3h2h4變步長二次偏導數(shù)（c）雙向插值法(i,j)(i+1,j)(i-1,j)01（2）第二類和第三類邊界條件對于點O過O點向邊界G做垂線PQ交邊界于Q，交網(wǎng)線段VR于P，OP=ah，PR=bh，VP=ch因為P一般不是節(jié)點，其值應當以點和P、R點的插值給出代入第二、三類邊界條件（2）第二類和第三類邊界條件對于點O過O點向邊界G做垂線圖中O與R重合圖中V與R點重合（2）第二類和第三類邊界條件圖中O與R重合圖中V與R點重合（2）第二類和第三類邊界條件2.4差分方程對于具體地球物理問題的偏微分方程組，利用上述差分格式，可以給出偏導數(shù)的微商近似，進一步得到差分方程組。設f(P)是內(nèi)部區(qū)域DI上定義的一個函數(shù)，設L(u)是一個微分算子，則以下表示了未知量u(P)的偏微分方程：邊界條件表示為以下方程：設和分別表示區(qū)域D的北部節(jié)點和邊界節(jié)點，則下式表示了以上偏微分方程的有限差分方程(finite-differenceequations,orfinite-differencescheme)：其中：以上差分方程是偏微分方程的有限差分近似，U是u的有限差分近似。差分方程要求U在所有節(jié)點上是u的很好的近似，并且方程所給出的有限差分近似解U是唯一的。2.4差分方程對于具體地球物理問題的偏微分方程組，利用上述例子：對流方程（雙曲型）的初值問題差分方程(13)(14)假定以上問題的解u(x,t)是充分光滑的，由Taylor級數(shù)展開有：例子：對流方程（雙曲型）的初值問題差分方程(13)(14)假利用（15）和（17）式，得到：如果u(x,t)是滿足方程（13）的光滑解，則代入（20），可以看出，偏微分方程（13）在(xj,tn)處可以近似地用下面的方程來代替：其中是u(xj,tn)的近似值。（21）式稱為逼近方程（13）的有限差分方程，簡稱差分方程。用到的節(jié)點如圖所示，（21）式可以改寫為便于計算的形式(20)(21)(22)利用（15）和（17）式，得到：如果u(x,t)是滿足方程（（22）式加上初始條件（14）的離散形式(22)（23）就可以按時間逐層推進，算出各層的值。這里的“層”表示在直線t=nτ上網(wǎng)格點的整體。差分方程（22）和初始條件的離散形式（23）結合在一起構成了一個差分格式。（22）給出了根據(jù)初始條件（23）來確定(j=0,±1,…)的一個算法，因此有時也稱差分方程（22）為一個差分格式。差分格式包含了初始條件、邊界條件的離散。由第n個時間層推進到第n+1個時間層時，（22）提供了逐點直接計算n+1時的表達式，因此稱（22）為顯式格式，在計算第n+1層時只用到了n層的數(shù)據(jù)，前后僅聯(lián)系到連個時間層，因此（21）式為兩層格式，更明確地稱其為兩層顯式格式。用（15）和（19），可以得到逼近方程（13）的另一差分格式其中λ為網(wǎng)格比。此格式也是兩層格式，稱為中心差分格式。（22）式加上初始條件（14）的離散形式(22)（23）就可用Taylor級數(shù)展開給出差分形式，用相應的差分形式逼近批未能微分方程（組），可以得到不同的差分格式，這一過程等價于用差商來近似微商得到相應的差分格式。不同的差分格式的精度和誤差不同。用Taylor級數(shù)展開給出差分形式，用相應的差分形式逼近批未例子：牛頓冷卻定律：溫度高于周圍環(huán)境的物體向周圍媒質(zhì)傳遞熱量逐漸冷卻時所遵循的規(guī)律。當物體表面與周圍存在溫度差時，單位時間從單位面積散失的熱量與溫度差成正比。Tair一階常微分方程的數(shù)值解如果不能簡單地通過積分求解，則需要利用數(shù)值方法求解。首先對時間和溫度進行離散：利用向前差分形式：得到以下的顯式差分格式：Tcap例子：牛頓冷卻定律：溫度高于周圍環(huán)境的物體向周圍媒質(zhì)傳遞熱量利用該差分格式，我們可以計算任一時間和函數(shù)f的溫度，但是隨著計算的進行，誤差O(dt)將會不斷積累。對于例子中的牛頓冷卻問題，我們想知道液體的溫度怎樣隨著時間和它與空氣的溫度差變化。設溫度差為τ為冷卻的時間尺度，這時常微分方程為初始條件：以上的差分格式為：該方程的解析解為：該有限差分格式的近似程度如何？怎樣選擇Δt才能得到穩(wěn)定解？利用該差分格式，我們可以計算任一時間和函數(shù)f的溫度，但是隨著差分方程該有限差分格式是否在Δt→0時收斂？它是否和解析解得性質(zhì)一樣？考慮初始條件得到因此，對于第j個時刻設t是到j時刻的總時間，則：對上式利用二項式展開其中：我們希望知道當dt→0時差分格式的結果，這時相當于j→∞差分方程該有限差分格式是否在Δt→0時收斂？它是否和解析解得因此：代入差分格式：上式就是解析解的Taylor展開。結論：對于牛頓冷卻問題，當時間步長Δt趨于零時，差分格式給出的數(shù)值解收斂于解析給出的嚴格解。解析解是單調(diào)減小函數(shù)，數(shù)值解的性質(zhì)怎樣？保證數(shù)值解單調(diào)減?。═j+1<Tj）需要什么條件？得到因此：代入差分格式：上式就是解析解的Taylor展開。結論：保證Tj+1<Tj的條件是：因此，數(shù)值解只對一定取值范圍的dt是單調(diào)減小的。如果數(shù)值解振蕩，但是結果會收斂如果數(shù)值解振蕩并且發(fā)散。因此，數(shù)值結果是有條件穩(wěn)定的。保證Tj+1<Tj的條件是：因此，數(shù)值解只對一定取值范圍的d差分方程fori=1:nt+1xt(i)=(i-1)*dt;T1(i)=t0*exp(-xt(i)/tau);endplot(xt,T1,'r.-');holdon

set(gca,'DataAspectRatio',[(max(xt)-min(xt))/(max(T)-min(T))/311]);xlabel('Time(s)','Fontname','timesnewroman','FontSize',14);

ylabel('Temperature','Fontname','timesnewroman','FontSize',14);%Malab-1Dclear;clc;figure('color','w');nt=8;%totaltimestepst0=1;%initialtemperaturetau=0.7;%timeconstantdt=1.25;%timeintervalT(1)=t0;fori=1:nt;xt(i)=(i-1)*dt;T(i+1)=T(i)-dt/tau*T(i);endxt(nt+1)=nt*dt;plot(xt,T,'b.-');holdon解析解差分解差分方程fori=1:nt+1%Malab-1D解析解差結果振蕩但是收斂，前期結果不準確結果收斂，前期結果不準確結果振蕩但是收斂，結果收斂，結果收斂，誤差逐漸縮小結果收斂，結果基本吻合解析解結果收斂，結果收斂，結果振蕩，不收斂，計算結果錯誤結果振蕩，不收斂，計算結果錯誤差分格式的性質(zhì)分析差分格式的性質(zhì)分析有限差分法基礎課件有限差分法基礎課件差分格式的性質(zhì)分析差分格式的性質(zhì)分析差分格式的性質(zhì)分析穩(wěn)定性（stability）：如果偏微分方程的嚴格解析解有界，差分格式給出的解也有界，稱該差分格式是穩(wěn)定的；如果差分格式給出的解是無界的，則稱該差分格式是不穩(wěn)定的。如果差分格式給出的解對于所有的時間步長和空間步長取值都是有界的，則稱該差分格式是無條件穩(wěn)定的；如果只是對時間步長和空間步長的部分取值有界，稱它是有條件穩(wěn)定的；如果對于所有的時間步長和空間步長取值都是無界的，則稱差分格式是無條件非穩(wěn)定的。穩(wěn)定性反映了差分格式在計算中控制誤差傳遞的能力，對偏微分方程有限差分方法研究具有重要意義。例子：對流方程（雙曲型）的初值問題差分格式的性質(zhì)分析穩(wěn)定性（stability）：如果偏微分方有限差分法基礎課件收斂性（convergence）：如果當時間和空間步長趨于零時，F(xiàn)DE解趨于PDE解，稱該差分格式是收斂的。如果則稱該差分格式是收斂的。對連續(xù)形式的偏微分方程進行有限差分離散時，差分格式對最終計算結果有重要影響。收斂性表示當差分網(wǎng)格無限細化時，差分方程的解是否具有無限逼近偏微分方程的解的能力。收斂性的判斷決定了一個差分格式能否被用來離散偏微分方程，不收斂的差分格式無疑是無法得到偏微分方程的真實解的。收斂性（convergence）：如果當時間和空間步長趨于零利用向前差分，得到以下差分格式：當時間步長dt趨近于零時，差分格式的近似解趨近于方程的解析解說明：該差分格式是收斂的。收斂性實例利用向前差分，得到以下差分格式：當時間步長dt趨近于零時，差很重要的一點：相容性是差分格式的性質(zhì)，它將差分格式和原始的偏微分方程聯(lián)系在一起。穩(wěn)定性和收斂性則是差分格式給出的數(shù)值解的性質(zhì)。一般來講，分析穩(wěn)定性和相容性比較容易，證明收斂性有時是很困難的數(shù)學問題。Lax等價定理（Laxequivalencetheorem）：如果逼近一個給定問題的差分格式是相容的，那么該差分格式的收斂性與穩(wěn)定性互為充分必要條件。該定理表明，對于一個具有相容性的差分格式，它的收斂性與穩(wěn)定性是等價的。如果格式不穩(wěn)定，則也不收斂。由于判斷差分格式的收斂性相對比較困難，該定理提供了通過判斷差分格式穩(wěn)定性來確定收斂性的方法。因此，一般不再討論收斂性問題，而是討論差分格式的穩(wěn)定性。該定理將收斂性與相容性和穩(wěn)定性聯(lián)系在一起，是很有用的一個定理。差分格式的性質(zhì)分析很重要的一點：相容性是差分格式的性質(zhì)，它將差分格式和原始的偏拉普拉斯方程U=0U=0U=100U=0拉普拉斯方程的有限差分法問題：2.5應用實例拉普拉斯方程U=0U=0U=100U=0拉普拉斯方程的問題：U=0U=0U=100U=0214365871091211141315U=0U=0U=100U=02143658710912111組建A和B矩陣，求解線性方程組得到X組建A和B矩陣，求解線性方程組得到X%Matlab2Dclear;clc;figure('color','w');

a=zeros(135,135);fori=1:135a(i,i)=1;end;fori=1:7a(15*i+1,15*i+2)=-0.25;a(15*i+1,15*i+16)=-0.25;a(15*i+1,15*i-14)=-0.25;endfori=1:7a(15*i+15,15*i+14)=-0.25;a(15*i+15,15*i+30)=-0.25;a(15*i+15,15*i)=-0.25;Enda(1,2)=-0.25;a(1,16)=-0.25;a(121,122)=-0.25;a(121,106)=-0.25;a(135,134)=-0.25;a(135,120)=-0.25;a(15,14)=-0.25;a(15,30)=-0.25;fori=2:14a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i+15)=-0.25;endfori=122:134a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i-15)=-0.25;endfori=1:7forj=2:14;a(15*i+j,15*i+j-1)=-0.25;a(15*i+j,15*i+j+1)=-0.25;a(15*i+j,15*i+j+15)=-0.25;a(15*i+j,15*i+j-15)=-0.25;endendb=a^(-1);