基于malab的偏微分方程箱在求解數(shù)學(xué)物理方程中的應(yīng)用

基于malab的偏微分方程箱在求解數(shù)學(xué)物理方程中的應(yīng)用

nb是世界上最流行的科技應(yīng)用之一。這是1984年由美國(guó)mathrs開發(fā)的一個(gè)軟件包理軟件。但是Matlab軟件在求解數(shù)學(xué)物理方程中的應(yīng)用并沒有受人們的足夠重視。在工程領(lǐng)域,一般都利用相關(guān)的其他應(yīng)用軟件(如ANSYS軟件等)或用VB、VC++等編程語(yǔ)言針對(duì)具體的數(shù)學(xué)物理方程編程,這樣需要消耗大量的人力、物力?；谝陨系默F(xiàn)實(shí)情況,本文主要討論MATLAB7.0的偏微分方程(PDE)工具箱在求解數(shù)學(xué)物理方程中的一些應(yīng)用。通過作者應(yīng)用研究經(jīng)驗(yàn)表明:MATLAB7.0的PDEToolbox可求解《數(shù)學(xué)物理方程》等的教材中列出來(lái)的大多數(shù)定解問題,只是這些教材中求出的是解析解,對(duì)于初邊值條件的可以用一個(gè)函數(shù)符號(hào)代替,而在用(MATLAB軟件求解時(shí)必須是已知的表達(dá)式或具體的數(shù)值。目前的MATLAB7.0的PDEToolbox提供了研究和求解空間二維偏微分方程問題的一個(gè)強(qiáng)大而又靈活實(shí)用的平臺(tái)。該工具箱能求解基本方程有橢圓型方程、拋物型方程、雙曲型方程、特征值方程、橢圓型方程組以及非線性橢圓型方程。1任意的工具函數(shù)Matlab工具箱就是一些M文件的集合,用戶可以修改工具箱中的函數(shù),更為重要的是用戶可以通過編制M文件來(lái)任意地添加工具箱中原來(lái)沒有的工具函數(shù)。許多的專業(yè)領(lǐng)域在MATLAB中都有自己的工具箱。下面僅就用Matlab7.0PDEToolBox求解數(shù)學(xué)物理方程需要注意的問題作簡(jiǎn)單介紹。1.1b接觸條件混合邊界條件橢圓型方程:-?·(c?u)+au=f(1)拋物型方程:d?u?t???(c?u)+au=f(2)d?u?t-??(c?u)+au=f(2)雙曲型方程:d?2u?t2???(c?u)+au=f(3)d?2u?t2-??(c?u)+au=f(3)特征值方程:-?·(c?u)+au=λdu(4)其中d是定義在Ω上的復(fù)變函數(shù),λ是待求的特征值。橢圓型方程組:??????????????(c11?u1)???(c12?u2)?????(c1n?un)+a11u1+a12u2+?+a1nun=f1???(c21?u1)???(c22?u2)?????(c2n?un)+a21u1+a22u2+?+a2nun=f2?????(cm1?u1)???(cm2?u2)?????(cmn?un)+am1u1+am2u2+?+amnun=fm(5){-??(c11?u1)-??(c12?u2)-?-??(c1n?un)+a11u1+a12u2+?+a1nun=f1-??(c21?u1)-??(c22?u2)-?-??(c2n?un)+a21u1+a22u2+?+a2nun=f2??-??(cm1?u1)-??(cm2?u2)-?-??(cmn?un)+am1u1+am2u2+?+amnun=fm(5)非線性橢圓方程-?·(c(u)?u)+a(u)u=f(6)邊界條件:①Dirichlet條件:hu=r(7)②Neumann條件:n·(c?u)+qu=g(8)要特別說明的是:對(duì)于特征值問題僅限于齊次條件;對(duì)于非線性情形,系數(shù)g,q,h和r可以依賴于u;對(duì)于拋物型方程和雙曲型方程,系數(shù)可以依賴于時(shí)間t。③Robin條件:就是未知函數(shù)在邊界上的法向微商與其本身的線性組合的值。④混合邊界條件:指既有Dirichlet條件又有Neumann條件,甚至還有的Robin條件的情形。初值條件:一般是含有時(shí)間(t)變量的拋物型方程和雙曲型方程等才涉及初值條件。在使用MATLAB7.0的PDEToolbox之前,要把原方程(組)先轉(zhuǎn)化成工具箱能識(shí)別的形式。1.2求解表的一般過程根據(jù)欲求解的方程(組)及相關(guān)的初邊值條件,用MATLAB7.0的PDEToolbox和使用函數(shù)命令及相關(guān)編程手段來(lái)求解時(shí)皆可總結(jié)為以下幾步:①根據(jù)未知變量的取值范圍畫出定解區(qū)域;②設(shè)置相關(guān)的初邊值條件;③設(shè)置方程類型;④網(wǎng)格剖分;⑤解方程。2使用示例下面僅用MATLAB7.0的PDEToolbox求解拋物型方程、雙曲型方程兩類典型的方程實(shí)例。2.1求解線性系統(tǒng)的解析解設(shè)有半徑為R的薄圓盤,其上下兩底絕熱,若圓盤邊界上的溫度保持為零度,且初始溫度為已知,求圓盤內(nèi)的瞬時(shí)溫度分布規(guī)律。這個(gè)問題可以歸結(jié)為求解下述的定解問題:???????a2?u?t=?2u?x2+?u?y2,x2+y2=R2,t>0ut=0=φ(x,y)ux2+y2=R2=0(9){a2?u?t=?2u?x2+?u?y2,x2+y2=R2,t>0ut=0=φ(x,y)ux2+y2=R2=0(9)這里假設(shè)R=0.5,a=1,ut=0=0.5,t>0。用MATLAB7.0的PDEToolbox經(jīng)過69步,140次函數(shù)調(diào)用,18次LU分解,139次的線性系統(tǒng)的求解求得本問題的色彩解如下:由圖1可以很直觀地看出該時(shí)刻圓盤的熱傳導(dǎo)規(guī)律:從圓盤的中央向四周逐漸傳熱,使得四周的溫度逐漸升高。由初邊界條件及相關(guān)的熱傳遞規(guī)律知其解的是符合實(shí)際情況的。也說明了其解的正確性。用分離變量法求解式(9)的解析解時(shí)可以引出著名Bessel方程和Helmholtz方程,其相應(yīng)的解析解是通過分離變量法轉(zhuǎn)化成求解Bessel方程的特征值問題來(lái)求解。詳細(xì)過程見參考文獻(xiàn)。2.2ca多通道波問題的解析解假設(shè)考慮如下的具體方程如下:??????????????????????2u?t2??2u=0,(x,y)∈{(x,y)?1<x<1,?1<y<1}ux=±0=0?u?ny=±1=0ut=0=arctan(cos(π2x))?u?tt=0=3sin(πx)esin(π2y)(10){?2u?t2-?2u=0,(x,y)∈{(x,y)-1<x<1,-1<y<1}ux=±0=0?u?ny=±1=0ut=0=arctan(cos(π2x))?u?tt=0=3sin(πx)esin(π2y)(10)用MATLAB7.0的PDEToolbox經(jīng)過106步,6次failedattempts,226函數(shù)調(diào)用,27次LU分解,最終求解了225次線性系統(tǒng)得本問題的色彩解如圖2所示。由圖2(b)可以直觀地看出在二維平面中波的傳播現(xiàn)象,如波峰、波谷等,如果區(qū)域選擇大一些波的轉(zhuǎn)播現(xiàn)象更能看得清楚。關(guān)于式(10)的解析解的問題,目前只能用降維法求解出相應(yīng)的Cauchy問題的解析解。由于篇幅限制和數(shù)值解本身是一個(gè)很大的數(shù)值矩陣的原因,兩個(gè)實(shí)例的相應(yīng)的數(shù)值解未在本文中列出。3數(shù)學(xué)物理方程的新方法、新出路與其他求解方法相比,此方法有計(jì)算速度快,計(jì)算方便,求解穩(wěn)定可靠、節(jié)省財(cái)力人力等優(yōu)點(diǎn),為廣大應(yīng)用科學(xué)的