




版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進(jìn)行舉報或認(rèn)領(lǐng)
文檔簡介
1、Possion 方程的差分方法課程名稱: 微分方程數(shù)值解 學(xué)生姓名: 張弘 一、問題描述二、問題分析I偏微分方程的數(shù)值解法主要有有限差分法和Galerkin有限元法。用差分法和有限元法將連續(xù)問題離散化的步驟是:1、對求解區(qū)域做網(wǎng)格剖分用有限個網(wǎng)格節(jié)點(diǎn)代替連續(xù)區(qū)域。2、微分算子離散化。3、把微分方程的定解問題化為線性代數(shù)方程組的求解問題。差分法和有限元方法的主要區(qū)別是離散化的第二步,差分法從定解問題的微分或積分形式出發(fā),用數(shù)值微商或數(shù)值積分公式到處相應(yīng)的線性代數(shù)方程組,有限元法從定解問題的變分形式出發(fā),用Ritz-Galerkin法導(dǎo)出相應(yīng)的線性代數(shù)方程組。差分法的基本問題有:(1) 對求解區(qū)域
2、做網(wǎng)格剖分一維情形為把區(qū)間分為等距或不等距的區(qū)間,二維情形是把區(qū)域分割成均勻或不均勻的矩形,其邊與坐標(biāo)軸平行,也可分割成小三角形或凸四邊形。(2) 構(gòu)造逼近微分方程定解問題的差分格式有兩種構(gòu)造差分格式的方法:直接差分化法和有限體積法。(3) 差分解的存在唯一性,收斂性和穩(wěn)定性研究(4) 差分方程的解法:逐次超松弛法,交替方向法,共軛梯度法。II(1)由題目可知,本題需要考慮矩形網(wǎng)的五點(diǎn)差分格式。題目給出的為第一邊值條件,將十字形圖形的中心放于坐標(biāo)原點(diǎn)處,如下圖所示:OS1G1S2由圖形可知,所給區(qū)域可以看出是由8個梯形G1通過旋轉(zhuǎn)、翻轉(zhuǎn)拼接而成。因此為了方便計算、減少計算量,只針對G1進(jìn)行網(wǎng)格
3、剖分,用5點(diǎn)差分格式進(jìn)行求解。但是由于G1是直角梯形,進(jìn)行網(wǎng)格剖分時會出現(xiàn)x軸方向網(wǎng)格點(diǎn)個數(shù)不同的現(xiàn)象,不利于有差分形成的系數(shù)矩陣的生成,所以將三角形S1和梯形G1合在一起形成一個矩形,其區(qū)域?yàn)?,3/2×0,1/2。如果采用等距差分,并且有hx=hy=h,設(shè)步長為h,x=0:h:3/2;y=0:h:1/2;nx=length(x)-1;%為所求區(qū)域中x軸上內(nèi)點(diǎn)的個數(shù)ny=length(y)-1; %為所求區(qū)域中y軸上內(nèi)點(diǎn)的個數(shù)對于原來的梯形G1來說,有網(wǎng)格內(nèi)點(diǎn)nx*ny-(ny2-my)/2對于矩形區(qū)域S1+G1而言,網(wǎng)格內(nèi)點(diǎn)數(shù)為nx*ny,所以在后面的程序中要比在梯形G1中多計算
4、了(ny2-my)/2 個點(diǎn)的函數(shù)值,但對程序效率的影響并不是很大,可以接受。以下具體說明只需在G1上求解poisson方程的原因所求方程為:設(shè)直線l是經(jīng)過原點(diǎn)o的任意一條直線,其方程為:y=kx設(shè)p(x,y)是區(qū)域內(nèi)任一點(diǎn),則其關(guān)于l對稱的點(diǎn)為p'(s,t)S=(k-1)2+|2y)/(k2+1) t=(2kx+(k2-1)y)/(k2+1)則同理可得:將u(s,t)代替u(x,y)得:Uxx+uyy=uss+utt=1其依然滿足上述方程。令=arctan(k)點(diǎn)p的橫坐標(biāo)x=rcos() y=rsin()則p關(guān)于直線l的對稱點(diǎn)為p'(rcos(2-),rsin(2-)由上述
5、證明可知u(p)=u(p')。由和p點(diǎn)的任意性知,對于函數(shù)u圖像上的任意一點(diǎn)p,其關(guān)于任意一條經(jīng)過原點(diǎn)直線l的對稱點(diǎn)p'都在u的圖像上,即u(+)=u(),即u關(guān)于原點(diǎn)是旋轉(zhuǎn)對稱的。當(dāng)l為x軸時,有u(x,y)=u(x,-y)l為y軸時,u(x,y)=u(-x,y)坐標(biāo)軸旋轉(zhuǎn)不改變方程的形式,所以函數(shù)在直角坐標(biāo)系中u關(guān)于原點(diǎn)是旋轉(zhuǎn)對稱的,又所求十字形區(qū)域關(guān)于x,y軸是軸對稱和關(guān)于原點(diǎn)中心對稱的,因此可通過直求解區(qū)域G1,就可以知道函數(shù)在整個十字形區(qū)域的圖像。(2)對S1+G1形成的矩形進(jìn)行正方形網(wǎng)格剖分,則求解Poisson方程的差分格式和化為如下形式:對于正則內(nèi)點(diǎn)其差分如下:
6、-uij=1/h2*(-u(i,j+1)-u(i-1,j)+4u(i,j)-u(i+1,j)-u(i,j-1)=fij(1)對于矩形區(qū)域S1+G1,nx=(xb-xa)/h-1ny=(yb-ya)/h-1按從左到右,從下到上的次序排列網(wǎng)點(diǎn)(ij):j=1,1inx;j=2,1inx;,;j=ny,1inx,定義向量Uh=u11,u21,unx-1;u1,nx-1,u2,nx-1,uny-1,nx-1T利用邊界條件可以將(1)式寫成如下形式:其中A=為nx*ny階矩陣,I為nx階單位矩陣,B為nx階單位矩陣。B=右端向量g的元素,依賴于邊值a和右端項(xiàng)f,顯然A是對稱正定矩陣,也是稀疏矩陣因此可以
7、采用逐次超松弛法,共軛梯度法和交替方向法萊求解,但為了方便程序設(shè)計采用了matlab的'運(yùn)算符來求解u。針對本題的Poisson方程,對S1+G1形成的矩形網(wǎng)格的五點(diǎn)差分的具體分析。對S1+G1形成的矩形五點(diǎn)差分的具體分析。ny+1ny3211 2 3 4 i nx-1 nx nx+1G1S2S11/2O3/21/2紅色線條圍成的區(qū)域?yàn)镚1,L為紅色斜邊,S1為L左側(cè)的三角形,S2為L右側(cè)的三角形。由對稱性知,S1和S2中關(guān)于L相互對稱的網(wǎng)格點(diǎn)其上U的函數(shù)值是相同的。按從左到右,從下到上的次序排列網(wǎng)點(diǎn)(ij)。(1)對于三角形S1中的網(wǎng)點(diǎn)有u(i,j), nyi>j1,有u(i,
8、j)-u(j,i)=0 所以S1中網(wǎng)點(diǎn)的差分就為:u(i,j)-u(j,i)=0(2)由于原點(diǎn)o點(diǎn)的特殊性,其鄰點(diǎn)中u(1,2)=u(1,-1)=u(-1,1)=u(2,1)所以其差分為:u(1,1)-4u(1,2)= h2*f(i,j)程序中語句:A(1:nx,nx+1:2*nx)=2*I; 和A(1,2)=-2;保障了上面差分方程的系數(shù)。(3)對于網(wǎng)格上最下面一行上除了原點(diǎn)o的所有正則網(wǎng)格內(nèi)點(diǎn),由對稱性得u(1,i)的鄰點(diǎn)中u(1,i-1)=u(1,i+1)所以其差分為:4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)=h2*f(i,j)對于右下角的非正則內(nèi)點(diǎn)u(1
9、,nx)其差分為:4u(i,j)-u(i-1,j)-2*u(i,j+1)=h2*f(i,j)程序中的相關(guān)語句為:A(1:nx,nx+1:2*nx)=2*I; A(nx+1:2*nx,nx+1:2*nx)=B;(4)對于G1中的所有正則內(nèi)點(diǎn),其差分為:4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j-1)-u(i,j+1)=h2*f(i,j)程序中相關(guān)語句:A(i*nx+1:(i+1)*nx,i*nx+1:(i+1)*nx)=B;A(i-1)*nx+1:i*nx,i*nx+1:(i+1)*nx)=I;A(i*nx+1:(i+1)*nx,(i-1)*nx+1:i*nx)=I;(5)對
10、于網(wǎng)格最后一列的所有非正則內(nèi)點(diǎn)u(:,nx),其差分為:4u(nx,j)-u(nx,j-1)-u(nx,j+1)-u(nx-1,j)=h2*f(i,j)(6)在求出了所有的內(nèi)點(diǎn)后,對于剩下的邊界點(diǎn)賦值有:U(ny+1,ny+1:nx)=0%最上面一行上的邊界點(diǎn)U(1:ny+1,nx+1)=0%最右側(cè)一列的邊界點(diǎn)u(ny+1,1:ny)=u(1:ny,ny+1);%三角形S1和S2邊界上的網(wǎng)點(diǎn),它們是S1的邊界點(diǎn),但是整個求解區(qū)域的內(nèi)點(diǎn)。三、程序設(shè)計及分析function poisson5spot(h)if nargin<1 %默認(rèn)步長h=0.02 h=0.02;endx=0:h:3/2;
11、y=0:h:1/2;nx=length(x)-1;%取區(qū)域的內(nèi)點(diǎn)ny=length(y)-1;%=構(gòu)造矩陣B=B=eye(nx,nx)*4;for i=1:nx-1 B(i,i+1)=-1;endfor i=2:nx B(i,i-1)=-1;endI=-eye(nx,nx);%=構(gòu)造系數(shù)矩陣A=A=zeros(nx*ny,nx*ny);A(1:nx,1:nx)=B;%由區(qū)域的對稱性,正方形網(wǎng)格最下面一行的差分形式%改為4u(i,j)-u(i-1,j)-u(i+1,j)-2*u(i,j+1)%因?yàn)閡(i,j+1)=u(i,j-1)A(1:nx,nx+1:2*nx)=2*I; A(nx+1:2*n
12、x,nx+1:2*nx)=B;%矩形網(wǎng)格左下角第一個點(diǎn)的差分形式改為:%4u(i,j)-u(i-1,j)-u(i+1,j)-u(i,j+1)-u(i,j-1)%=4u(i,j)-2u(i,j+1)-2u(i+1,j)A(1,2)=-2;%為了方便,本程序往梯形中增加了一個三角形,以方便編程A(nx+1:2*nx,1:nx)=I;for i=2:ny-1 A(i*nx+1:(i+1)*nx,i*nx+1:(i+1)*nx)=B; A(i-1)*nx+1:i*nx,i*nx+1:(i+1)*nx)=I; A(i*nx+1:(i+1)*nx,(i-1)*nx+1:i*nx)=I;end%=%bi=h
13、*h*f;右端向量b=h*h*ones(nx*ny,1);%由于對稱性,左側(cè)三角形中有u(i,j)-u(j,i)=0 i>j%因此令A(yù)(i,j)=1 A(j,i)=-1%所以在本程序中多計算了(ny2-my)/2 個點(diǎn)的函數(shù)值%但對程序的影響并不是很大for i=2:ny A(i-1)*nx+1:(i-1)*nx+i-1,:)=0; for j=1:i-1 A(i-1)*nx+j,(i-1)*nx+j)=1; A(i-1)*nx+j,(j-1)*nx+i)=-1; b(i-1)*nx+j)=0; endendx=Ab;%為了方便采用左除求解網(wǎng)格點(diǎn)上的函數(shù)值%x=gmres(A,b,100
14、);%按順序?qū)賦值給uu=zeros(ny+1,nx+1);for i=1:ny for j=1:nx u(i,j)=x(i-1)*nx+j); endend%根據(jù)對稱性,給網(wǎng)格最上面一行的點(diǎn)賦值u(ny+1,1:ny)=u(1:ny,ny+1); %=作Poisson方程在區(qū)域上的圖形=x,y=meshgrid(0:h:3/2,0:h:1/2);hold onsurf(x,y,u)%11第一象限的第一塊區(qū)域,下面的以此類推surf(y,x,u)%12surf(-y,x,u)%21surf(-x,y,u);%22surf(-x,-y,u)%31surf(-y,-x,u)%32surf(y,-x,u)%41surf(x,-y,u);%42 四、實(shí)驗(yàn)結(jié)果1在區(qū)域G1上求解后的圖形顯示:圖形表示了在S1+G1區(qū)域上的函數(shù)圖像,而不是單純的G1上的函數(shù)圖像。2通過拼接后圖形顯示如下:由上圖可知,除了邊界點(diǎn)外網(wǎng)格點(diǎn)上的函數(shù)值都有u(i,j)>0.Lhuij>0,對任意(xi,yj)Gh,uij不可能在內(nèi)點(diǎn)取得負(fù)的極小,與極值定理相符合。3、翻轉(zhuǎn)后圖形如下:4 網(wǎng)格間距h=0.01時得到
溫馨提示
- 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
- 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
- 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
- 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
- 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負(fù)責(zé)。
- 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請與我們聯(lián)系,我們立即糾正。
- 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。
最新文檔
- 無線基站維護(hù)培訓(xùn)課件
- 抖音商戶短視頻創(chuàng)意提案評審制度
- BWA-6047-生命科學(xué)試劑-MCE
- 江蘇省興化市顧莊區(qū)三校2024-2025學(xué)年七上數(shù)學(xué)期末教學(xué)質(zhì)量檢測模擬試題含解析
- 美發(fā)培訓(xùn)卷杠課件
- 國際多式聯(lián)運(yùn)操作規(guī)范與風(fēng)險管理
- 航空行業(yè)三年發(fā)展報告:國際與國內(nèi)市場的比較研究
- 2024-2025學(xué)年浙江省杭州市濱江區(qū)數(shù)學(xué)七年級第一學(xué)期期末調(diào)研試題含解析
- 云南司法警官職業(yè)學(xué)院《國畫山水》2023-2024學(xué)年第一學(xué)期期末試卷
- 河道垃圾清理管理辦法
- 2025屆黑龍江省大慶中學(xué)九上化學(xué)期末聯(lián)考試題含解析
- 20濕性愈合功能性敷料的種類與敷料選擇
- 燃?xì)鈭缶餍袠I(yè)發(fā)展分析及投資戰(zhàn)略研究報告2025-2028版
- 2025年中國扭蛋行業(yè)市場全景分析及前景機(jī)遇研判報告
- 2025至2030中國現(xiàn)金處理中心行業(yè)發(fā)展趨勢分析與未來投資戰(zhàn)略咨詢研究報告
- 小學(xué)音標(biāo)題目及答案
- 2024年宿州蕭縣縣直事業(yè)單位招聘真題
- 美好生活大調(diào)查:中國居民消費(fèi)特點(diǎn)及趨勢報告(2025年度)
- 2025河南省豫地科技集團(tuán)有限公司社會招聘169人筆試參考題庫附帶答案詳解
- 快遞分揀人力承包協(xié)議書
- Q-GDW10162-2025 輸電桿塔固定式防墜落裝置技術(shù)規(guī)范
評論
0/150
提交評論