數(shù)值模擬1(FDM)課件

1、地球物理數(shù)值模擬方法地球探測(cè)與信息技術(shù)專業(yè)研究生課程一、 課程性質(zhì)、目的和任務(wù)地球物理數(shù)值模擬是地球探測(cè)與信息技術(shù)專業(yè)應(yīng)用地球物理方向的一門(mén)必修課。通過(guò)學(xué)習(xí)，使學(xué)生了解和掌握用有限元、邊界元和有限差分法解地球物理中的邊值問(wèn)題。二、 教學(xué)要求、形式和方法要求學(xué)完以后，學(xué)生能夠用有限元、邊界元、有限差分中的一種方法編制某一實(shí)際問(wèn)題的數(shù)值模擬計(jì)算程序，并撰寫(xiě)讀書(shū)報(bào)告。在教學(xué)形式上以自學(xué)為主，課堂教學(xué)為輔，教學(xué)方法上則采用教師引導(dǎo)性講學(xué)與學(xué)生結(jié)合問(wèn)題的自學(xué)。 三、 教學(xué)內(nèi)容第四講：邊界元結(jié)合地形影響的邊界元數(shù)值模擬講述：積分方程的導(dǎo)出，單元剖分，插值方法，數(shù)值積分方法，線性方程解法。第五講：成績(jī)?cè)u(píng)定對(duì)

2、所學(xué)知識(shí)，特別是自己所研究解決的問(wèn)題作學(xué)術(shù)報(bào)告。 四、教材及參考書(shū)目1、地球物理中的有限單元法，徐世浙，科學(xué)出版社，1994；2、地球物理中的邊界單元法，徐世浙，科學(xué)出版社，1995；3、電子計(jì)算機(jī)在電法勘探中的應(yīng)用，羅延鐘，張桂青，武漢地質(zhì)學(xué)院出版社，1987；4、電法勘探數(shù)值模擬技術(shù)，周 襄，鐘本 等編著，四川科學(xué)出版社，1986；5、電法勘探用邊界元法，田憲謨，黃蘭珍，地質(zhì)出版社，1990；6、The Finite Element Method, O.C.Zienkiewicz, R.L.Taylor, Mcgraw-Hill International Editions, 1989；7

3、、阮百堯，熊彬，徐世浙，2001，三維地電斷面電阻率測(cè)深有限元數(shù)值模擬，地球科學(xué)，26（1），73-77；8、阮百堯，2001，視電阻率對(duì)模型電阻率的偏導(dǎo)數(shù)矩陣計(jì)算方法，地質(zhì)勘探，37（6），39-41；9、阮百堯，2001，三角單元剖分電導(dǎo)率分塊連續(xù)變化點(diǎn)源二維電場(chǎng)有限元數(shù)值模擬，廣西科學(xué)，8（1），1-3；10、阮百堯，徐世浙，1998，電導(dǎo)率分塊線性二維地電斷面電阻率測(cè)深有限元數(shù)值模擬，地球科學(xué),23,303-307。第一講：概論正演問(wèn)題：已知場(chǎng)源和地下物性分布，計(jì)算地球物理場(chǎng)的空間和時(shí)間分布。對(duì)直流電阻率法而言，正演問(wèn)題是已知供電電極位置、電流大小，地下介質(zhì)的電阻率（或電導(dǎo)率）分布，求

4、地表和地下的電位。第一講：概論正演的作用：勘探設(shè)計(jì)的依據(jù)：針對(duì)不同的地質(zhì)條件、地質(zhì)目標(biāo)，確定用何種地球物理方法、測(cè)量的方式、測(cè)網(wǎng)、極距的大小，了解可能得到的異常；實(shí)測(cè)數(shù)據(jù)解釋的基礎(chǔ)：定性分析、一維、二維和三維的解釋。剖面I剖面II剖面I還是剖面II？第一講：概論正演的方法物理模擬法：水槽、土槽、導(dǎo)電紙、電阻網(wǎng)絡(luò)數(shù)學(xué)模擬法：數(shù)值計(jì)算法：有限差分法、有限元法、積分方程法（邊界元法）解析法：分離變量法、鏡像法、類比法數(shù)學(xué)模擬法：大部分的地球物理過(guò)程都能以數(shù)學(xué)的形式加以描述。 給定模型參數(shù)，數(shù)據(jù)觀測(cè)方式，建立數(shù)學(xué)物理模型，如：積分方程；微分方程、初始條件、邊界條件變分方程用解析或數(shù)值模擬方法求解數(shù)學(xué)物

5、理模型，得到預(yù)測(cè)數(shù)據(jù)。數(shù)值模擬方法常用的數(shù)值模擬方法有：有限差分法（FiniteDifferenceMethodFDM）有限元法（Finite Element MethodFEM ）積分方程法 邊界元法（Boundary Element MethodBEM）有限差分方法(FDM)是計(jì)算機(jī)數(shù)值模擬最早采用的方法，至今仍被廣泛運(yùn)用。該方法將求解域劃分為差分網(wǎng)格，用有限個(gè)網(wǎng)格節(jié)點(diǎn)代替連續(xù)的求解域。有限差分法以Taylor級(jí)數(shù)展開(kāi)等方法，把控制方程中的導(dǎo)數(shù)用網(wǎng)格節(jié)點(diǎn)上的函數(shù)值的差商代替進(jìn)行離散，從而建立以網(wǎng)格節(jié)點(diǎn)上的值為未知數(shù)的代數(shù)方程組。該方法是一種直接將微分問(wèn)題變?yōu)榇鷶?shù)問(wèn)題的近似數(shù)值解法，數(shù)學(xué)概念

6、直觀，表達(dá)簡(jiǎn)單，是發(fā)展較早且比較成熟的數(shù)值方法。有限差分法的缺點(diǎn)是必需進(jìn)行整個(gè)區(qū)域的剖分，并且要求網(wǎng)格比較規(guī)則，空間網(wǎng)格最好為直角網(wǎng)格。有限差分法的特點(diǎn)有限元法的特點(diǎn)有限元方法（FEM）最早應(yīng)用于結(jié)構(gòu)力學(xué)，其基本求解思想是借助于變分原理或加權(quán)余量法，將微分方程轉(zhuǎn)變?yōu)榉汉瘶O值問(wèn)題，然后把計(jì)算域劃分為有限個(gè)互不重疊的單元，在每個(gè)單元內(nèi)，將連續(xù)函數(shù)的泛函，離散成單元節(jié)點(diǎn)上函數(shù)值的泛函。根據(jù)泛函取極值的條件，得到各節(jié)點(diǎn)的函數(shù)值應(yīng)滿足的線性代數(shù)方程組。解代數(shù)方程組，得到各節(jié)點(diǎn)的函數(shù)值。FEM法的優(yōu)點(diǎn)是求解規(guī)范（甚至可用軟件產(chǎn)生FEM源碼），并且可以模擬所有形體。FEM的缺點(diǎn)也是必需進(jìn)行整個(gè)區(qū)域的剖分，與

7、FDM相比，相對(duì)比較復(fù)雜。積分方程法（邊界元法）的特點(diǎn)邊界元法（BEM）通過(guò)格林公式，將偏微分方程表示的邊值問(wèn)題轉(zhuǎn)變成區(qū)域邊界上的積分方程，然后將邊界劃分為有限個(gè)互不重疊的單元，進(jìn)行數(shù)值積分。BEM是在有限單元法以后發(fā)展起來(lái)的一種數(shù)值計(jì)算方法，可以說(shuō)是解邊界積分方程的有限元法。BEM的優(yōu)點(diǎn)是將區(qū)域問(wèn)題轉(zhuǎn)變?yōu)閰^(qū)域邊界問(wèn)題，使求解問(wèn)題的維數(shù)降低，三維變?yōu)槎S、二維變?yōu)橐痪S。BEM的缺點(diǎn)是不規(guī)范，當(dāng)區(qū)域內(nèi)有不同數(shù)量的介質(zhì)體時(shí)，需重新編制程序；另外還有格林公式中的基本解很難找，線性方程組的系數(shù)矩陣為滿秩矩陣，這些都導(dǎo)致BEM使用不如其他兩種方法廣泛。數(shù)學(xué)補(bǔ)充：關(guān)于偏微分方程的一些基本概念幾個(gè)典型方程拉

8、普拉斯算符哈密頓算符幾個(gè)典型方程波動(dòng)方程在一些聲學(xué)、電磁學(xué)、光學(xué)和力學(xué)的波動(dòng)問(wèn)題中常常出現(xiàn)這類方程。擴(kuò)散方程、熱傳導(dǎo)方程其中， 是擴(kuò)散過(guò)程中某種物質(zhì)的濃度，或是固體的傳熱過(guò)程中在x處、t時(shí)刻的溫度。定解條件對(duì)于一個(gè)n階常微分方程，常?？梢詫⑵浣鈱?xiě)成依賴于n個(gè)任意常數(shù)的通解形式。但是對(duì)于偏微分方程，情況要復(fù)雜得多，一般很難用通解的形式表示。我們都是特定條件下求方程的解，這樣的條件稱為定解條件。如果一個(gè)偏微分方程定解問(wèn)題在某個(gè)函數(shù)集合中存在唯一的解，而且在定解條件的原始資料有微小變化時(shí)，在某種義下解也僅有微小的變化，我們說(shuō)解關(guān)于定解條件是穩(wěn)定的。如果一個(gè)定解問(wèn)題解的存在性、唯一性和穩(wěn)定性都成立，就

9、稱定解問(wèn)題是適定的。我們所討論的定解問(wèn)題都是適定的。第二講：有限差分法差分一維有限差分二維有限差分時(shí)間域有限差分例題第1章 差分一階差分xx+xx-x第2章 一維波動(dòng)方程的有限差分法把上面得到的差分算子代替彈性波動(dòng)方程中的微分算子，波動(dòng)方程可以近似表示為：或者這個(gè)有限差分方程的截?cái)嗾`差是初始條件為：邊界條件為：有限差分方程和離散形式的初值邊界條件一起構(gòu)成了波動(dòng)方程的差分格式。如果已知在時(shí)間 和 ，空間 和 處的位移值，則可以根據(jù)差分波動(dòng)方程中時(shí)間 和空間 處的位移值的顯式表達(dá)式，通過(guò)簡(jiǎn)單數(shù)學(xué)運(yùn)算實(shí)現(xiàn)對(duì)整個(gè)彈性波波場(chǎng)的時(shí)間空間逐層推進(jìn)計(jì)算。這種有限差分格式稱為顯式格式。在解題過(guò)程中，需要利用初始

10、條件，設(shè) 時(shí)， 用中心差分可表示為：于是得 把它代入上述方程中的最后一項(xiàng)，并取 j0， 即對(duì)應(yīng)t0 ，得出舉例初位移不為零，初速為零時(shí)的例子：設(shè)初位移為：其余令 ， ，將微分方程改寫(xiě)成差分方程，即有其中，初始條件可以表示成此處，上標(biāo)為一時(shí)為t0網(wǎng)格剖分有限差分方法用差分算子代替微分算子，將連續(xù)的波動(dòng)方程化解為離散形式。因此，連續(xù)的波動(dòng)方程求解域也要進(jìn)行離散化，通常是對(duì)求解域作網(wǎng)格剖分。圖中給出了一維波動(dòng)方程的時(shí)間和空間的網(wǎng)格剖分。規(guī)則網(wǎng)格和不規(guī)則網(wǎng)格通常在有限差分方法中，空間步長(zhǎng)一般取相等距離，這種剖分方法稱為規(guī)則網(wǎng)格或者均勻網(wǎng)格剖分。在二維彈性波動(dòng)方程有限差分離散中 ， 兩個(gè)空間坐標(biāo)方向的空

11、間步長(zhǎng)，通常也取相等的距離。如果 ， 分別取不同的值，則可以稱為矩形不規(guī)則網(wǎng)格。clearN=4010; dx=0.0024;dt=0.0005; c=dt*dt/dx/dx;x=linspace(0,1,420);u(1:420,1)=0;u(181:240,1)=0.05*sin(pi*x(181:240)*7);u(2:419,2)=u(2:419,1)+c/2*(u(3:420,1)-2*u(2:419,1)+u(1:418,1);h=plot(x,u(:,1),linewidth,3);axis(0,1,-0.05,0.05);set(h,EraseMode,xor,makersiz

12、e,18)for k=2:N set(h,Xdata,x,Ydata,u(:,2); drawnow; u(2:419,3)=2*U(2:419,2)-U(2:419,1)+C*(U(3:420,2). -2*U(2:419,2)+U(1:418,2); u(2:419,1)=u(2:419,2); u(2:419,2)=u(2:419,3);end其中，假設(shè)在t0前一時(shí)刻波場(chǎng)位移是 ，根據(jù)中心差分公式得到具有二階精度 的差分格式：在上面介紹的初始條件的離散公式中，截?cái)嗾`差是O(t)，誤差階數(shù)低于整個(gè)差分方程的截?cái)嗾`差 ，誤差精度降低，可以通過(guò)在時(shí)間軸上增加虛擬網(wǎng)格點(diǎn)的方法提高精度。 作業(yè)將課

13、中的例題看懂，然后將這個(gè)matlab程序改寫(xiě)成fortran程序。有限差分格式相容性、收斂性和穩(wěn)定性第2章 一維有限差分一維泊松方程及定解問(wèn)題可描述許多簡(jiǎn)單的物理現(xiàn)象 (例如)：彈性棒之變形張力作用下弦之變形棒之溫度分布第2章 一維有限差分一維泊松方程式解的性質(zhì)解 總是存在 總是較數(shù)據(jù) “平滑”若對(duì)于所有 , , 則對(duì)于所有 , 已知 則解 為唯一有限差分離散（網(wǎng)格剖分）第2章 一維有限差分將區(qū)間 細(xì)分成 個(gè)相等的子區(qū)間，則其中有限差分近似（代替導(dǎo)數(shù)）第2章 一維有限差分3、有限差分方程（代數(shù)方程）第2章 一維有限差分則 3、有限差分方程（代數(shù)方程）第2章 一維有限差分(對(duì)稱)第2章 一維有限

14、差分3、有限差分方程（代數(shù)方程） 可以證明A為對(duì)稱正定矩陣，線性方程解存在且唯一。對(duì)任意因此對(duì)任意( 為對(duì)稱正定)存在且唯一例題：第2章 一維有限差分其中令例題：第2章 一維有限差分有限差分收斂性？第2章 一維有限差分1. 離散解 是否仍保有連續(xù)解 之定性性質(zhì)？2. 當(dāng) 時(shí)解是否更精確？3.對(duì)于 是否可使任意變小？第2章 一維有限差分有限差分收斂性？離散穩(wěn)定性令可證：有限差分收斂性？舍位誤差第2章 一維有限差分對(duì)任意 可證明令第2章 一維有限差分有限差分收斂性？舍位誤差令 為離散誤差相減得且第2章 一維有限差分有限差分收斂性？舍位誤差第2章 一維有限差分利用離散穩(wěn)定性估計(jì)或有限差分收斂性？舍位

15、誤差第2章 一維有限差分例題：第2章 一維有限差分例題：漸近于，第3章 二維有限差分二維泊松方程及定解問(wèn)題在于第3章 二維有限差分有限差分離散（網(wǎng)格剖分）第3章 二維有限差分有限差分近似（代替導(dǎo)數(shù)）第3章 二維有限差分3、有限差分方程（代數(shù)方程）令 .例：第3章 二維有限差分第3章 二維有限差分例：具有帶狀結(jié)構(gòu)，頻寬：第3章 二維有限差分方程式編號(hào)為 之分量第3章 二維有限差分分塊矩陣分塊三對(duì)角線矩陣第3章 二維有限差分分塊定義第3章 二維有限差分3、有限差分方程（代數(shù)方程） 可以證明A為對(duì)稱正定矩陣，線性方程解存在且唯一。對(duì)于任何，則( 為 對(duì)稱正定): 存在且唯一還可證明：第3章 二維有限差分有限差分收斂性