二維熱傳導(dǎo)方程有限差分法的matlab實(shí)現(xiàn)(編輯修改稿)
2024-09-02 07:56 本頁面
【文章內(nèi)容簡(jiǎn)介】 式。用一維向后差分格式的直接推廣是 （）此格式的截?cái)嗾`差仍為，仍用Fourier方法來分析這個(gè)格式的穩(wěn)定性，仿前可以得出其增長(zhǎng)因子是：因此有，即差分格式（）是絕對(duì)穩(wěn)定的。為了提高精度，對(duì)微分方程（）也可以用CrankNicolson型差分格式，這也是一維問題的直接推廣。其格式可寫為 （）這也是二階精度格式，其增長(zhǎng)因子是因此，對(duì)任何都有，所以()式也是絕對(duì)穩(wěn)定的?，F(xiàn)在考慮一下隱式格式（）式和（）式的求解方法。我們知道，在一位格式形成的方程組是系數(shù)矩陣為三對(duì)角矩陣的線性代數(shù)方程組，因此用追趕法很容易求解。而對(duì)于（）式和（）式導(dǎo)出的系數(shù)矩陣不是三對(duì)角矩陣，因此求解就不容易了。我們對(duì)于顯式格式和隱式格式的分析知道，在實(shí)際使用上都受到限制，因此構(gòu)造每層計(jì)算量不大的絕對(duì)穩(wěn)定的格式就成為一個(gè)具有現(xiàn)實(shí)意義并很有興趣的問題。在一維中，隱式格式是絕對(duì)穩(wěn)定并可用追趕法很容易求解。由此產(chǎn)生了下面將使用的交替方程隱式格式。它具有絕對(duì)穩(wěn)定、容易求解和有相當(dāng)精度的特點(diǎn)。 建立方程組我們?cè)跇?gòu)造微分方程（）的隱式格式和顯式格式中，對(duì)和做了同樣的處理，即同時(shí)在第層或第層取值。為了構(gòu)造一維形式的隱式格式，對(duì)二階導(dǎo)數(shù)用在第層上用未知的二階中心差分來代替，而則用在第層上用已知的二階中心差分來代替，這樣得到的方程組在僅方向是隱式的。比較容易求解，用追趕法就可以了。同理，為對(duì)稱起見，在下一時(shí)間層上重復(fù)上述步驟，即又僅在方向是隱式的，對(duì)方向是顯式的。這樣相鄰的兩個(gè)時(shí)間層合并起來構(gòu)成一個(gè)差分格式。故用多次追趕法就可以解出了。我們解向后差分格式的方程。令，則（）式變?yōu)? （）（）在x方向是隱式時(shí)，代入（），變形為 （）（）化為矩陣形式： （）（）在y方向是隱式時(shí)，代入（），可變形為 （）（）化為矩陣形式： （）第5章 MATLAB編程 追趕法對(duì)于前面求出的（）和（）矩陣，我們稱之為三對(duì)角矩陣。在一些實(shí)際問題中，例如解常微分方程的邊值問題，接熱傳導(dǎo)方程一級(jí)船體數(shù)學(xué)放樣中建立三次樣條函數(shù)等，都會(huì)要求解系數(shù)矩陣為對(duì)角占優(yōu)的三對(duì)角方程組，解三對(duì)角矩陣我們有一種特殊的方法叫追趕法。例如：求解線性方程組，其中為三對(duì)角矩陣。首先對(duì)矩陣做如下三角分解：其中為下三角矩陣，為上三角矩陣，且主對(duì)角線元素全為1.求解等價(jià)于求解兩個(gè)三角形方程組。先求解的解，再求解的解。追趕法公式實(shí)際上就是把高斯消去法用到求解三對(duì)角方程組上去的結(jié)果。這是由于特別簡(jiǎn)單，因此使得求解的計(jì)算公式非常簡(jiǎn)單，而且計(jì)算量很小。另外追趕法的計(jì)算公式中不會(huì)出現(xiàn)中間結(jié)果數(shù)量級(jí)的巨大增長(zhǎng)和舍入誤差的嚴(yán)重積累。最后在MATLAB中編程實(shí)現(xiàn)。追趕法程序：見附錄。 用追趕法求解線性方程組 解：在MATLAB命令窗口中輸入求解程序： A=[2 1 0 0。1 3 2 0。0 2 4 3。0 0 3 5]A = 2 1 0 0 1 3 2 00 2 4 3 0 0 3 5 f=[6 1 2 1]’。x=chase(A,f)輸出計(jì)算結(jié)果為：x = 9 12 1 0.即追趕法求出了方程組的解為： 二維熱傳導(dǎo)方程有限差分解的MATLAB編程程序見附錄。第6章 數(shù)值舉例分析 算 例為了便于觀察和表示微分方程的解，我們將最終解用圖像表示出來。 求滿足下列條件的二維熱傳導(dǎo)方程的解：在MATLAB命令窗口中輸入求解程序：a=。u_xy0=inline(39。039。,39。x39。,39。y39。)。u_xyt=inline(39。x^2*sin(y)y^2*sin(x)39。,39。x39。,39。y39。,39。t39。)。D=[0,5,0,5]。T=1000。Mx=50。My=50。N=50。[u,x,y,t]=sjy(a,D,T,u_xy0,u_xyt,Mx,My,N)。mesh(x,y,u)xlabel(39。x39。)ylabel(39。y39。)zlabel(39。U39。)輸出計(jì)算結(jié)果為：rx = +003ry =+003此解即為差分方程的近似解。 求解如下二維熱傳導(dǎo)方程在MATLAB命令窗口中輸入求解程序：a=10。u_xy0=inline(39。039。,39。x39。,39。y39。)。u_xyt=inline(39。x*sin(pi*
環(huán)評(píng)公示相關(guān)推薦
鄂ICP備17016276號(hào)-1