微分方程數(shù)值解法編程作業(yè)三_第1頁(yè)
微分方程數(shù)值解法編程作業(yè)三_第2頁(yè)
微分方程數(shù)值解法編程作業(yè)三_第3頁(yè)
微分方程數(shù)值解法編程作業(yè)三_第4頁(yè)
微分方程數(shù)值解法編程作業(yè)三_第5頁(yè)
已閱讀5頁(yè)，還剩4頁(yè)未讀， 繼續(xù)免費(fèi)閱讀

下載本文檔

版權(quán)說(shuō)明：本文檔由用戶提供并上傳，收益歸屬內(nèi)容提供方，若內(nèi)容存在侵權(quán)，請(qǐng)進(jìn)行舉報(bào)或認(rèn)領(lǐng)

文檔簡(jiǎn)介

1、微分方程數(shù)值解法編程作業(yè)三用數(shù)值方法求解正方形域上的Poisson方程邊值問(wèn)題解: （一）用Gauss-seidel迭代法求解線性方程組Au=f，具體程序和結(jié)果如下：程序：function u,k,T=tygs(n)h=1/(n+1);f(2:n+1,2:n+1)=2*h2;u=zeros(n+2,n+2);e=0.000000001;tic;for k=1:1000 %迭代求解 er=0; for j=2:n+1 for i=2:n+1 Ub=u(i,j); u(i,j)=(u(i-1,j)+u(i+1,j)+u(i,j-1)+u(i,j+1)+f(i,j)/4;er=er+abs(Ub-u

2、(i,j); %估計(jì)誤差end endif er/n2<e,break;end %判斷是否達(dá)到計(jì)算精度，如果達(dá)到則退出循環(huán)endT=tock結(jié)果：tygs(8)T = 0.0028k = 131ans = 0 0 0 0 0 0 0 0 0 0 0 0.0300 0.0476 0.0577 0.0623 0.0623 0.0577 0.0476 0.0300 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0577 0.0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0

3、577 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0577 0.0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0577 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0300 0.0476 0.0577 0.0623 0.0623 0.0577 0.0476 0.0

4、300 0 0 0 0 0 0 0 0 0 0 0（二）用塊Gauss-Seidel迭代法求解線性方程組Au=f，具體程序及結(jié)果如下：程序：functionu,k,T=tygsk(n)% u: 方程組的解； k: 迭代次數(shù)； n: 非邊界點(diǎn)數(shù)% a: 方程組系數(shù)矩陣的下對(duì)角線元素； b: 方程組系數(shù)矩陣的主對(duì)角線元素；% c: 方程組系數(shù)矩陣的上對(duì)角線元素；d: 追趕法所求方程的右端向量；% e: 允許誤差界； er:迭代誤差；h=1/(n+1);f(2:n+1,2:n+1)=2*h2;a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(

5、n+2,n+2); e=0.000000001;tic;for k=1:1000 er=0; for j=2:n+1 Ub=u(:,j); d(1:n)=f(2:n+1,j)+u(2:n+1,j-1)+u(2:n+1,j+1) ; % x=zg(a,b,c,d); % 用追趕法求解u(2:n+1,j)=x' er=er+abs(Ub-u(:,j); end if er/n2<e ,break;end % 判斷是否達(dá)到計(jì)算精度，如果達(dá)到則退出循環(huán)endT=tocksurf(u)function x=zg(a,b,c,d)n=length(b); %LU分解。u(1)=b(1); f

6、or k=2:n if u(k-1)=0 ,D=0,return; end l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1);end% 追趕法求解之追過(guò)程，求解Ly=d。y(1)=d(1); for k=2:n y(k)=d(k)-l(k)*y(k-1);end if u(n)=0 ,D=0,return; endx(n)=y(n)/u(n); for k=n-1:-1:1 x(k)=(y(k)-c(k)*x(k+1)/u(k); end結(jié)果：tygsk(8)T = 0.1563k = 64ans = 0 0 0 0 0 0 0 0 0 0 0 0.0300 0.

7、0476 0.0577 0.0623 0.0623 0.0577 0.0476 0.0300 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0577 0.0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0577 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0577 0.

8、0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0577 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0300 0.0476 0.0577 0.0623 0.0623 0.0577 0.0476 0.0300 0 0 0 0 0 0 0 0 0 0 0（三）用Jacobi迭代法求解線性方程組Au=f，具體程序和結(jié)果如下：程序：function u,k,T=tyja(n)h=1/(n+1);f(2:n+1,2:n+1)=2*h2;u=zeros(n+2,n+2);e=0.0

9、00000001;tic;for k=1:1000 %迭代求解 er=0; t=u; for j=2:n+1 for i=2:n+1 Ub=u(i,j); u(i,j)=(t(i-1,j)+u(i+1,j)+t(i,j-1)+u(i,j+1)+f(i,j)/4;er=er+abs(Ub-u(i,j); %估計(jì)誤差end endif er/n2<e,break;end %判斷是否達(dá)到計(jì)算精度，如果達(dá)到則退出循環(huán)endT=tock結(jié)果：tyja(8)T = 0.0058k = 249ans = 0 0 0 0 0 0 0 0 0 0 0 0.0300 0.0476 0.0577 0.0623

10、 0.0623 0.0577 0.0476 0.0300 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0577 0.0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0577 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0577 0.0961 0.1195 0.1305

11、 0.1305 0.1195 0.0961 0.0577 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0300 0.0476 0.0577 0.0623 0.0623 0.0577 0.0476 0.0300 0 0 0 0 0 0 0 0 0 0 0（四）用塊Jacobi迭代法求解線性方程組Au=f，具體程序及結(jié)果如下：程序：functionu,k,T=tyjak(n)% u: 方程組的解； k: 迭代次數(shù)； n: 非邊界點(diǎn)數(shù)% a: 方程組系數(shù)矩陣的下對(duì)角線元素； b: 方程組系數(shù)矩陣的主對(duì)角線元素；%

12、 c: 方程組系數(shù)矩陣的上對(duì)角線元素；d: 追趕法所求方程的右端向量；% e: 允許誤差界； er:迭代誤差；h=1/(n+1);f(2:n+1,2:n+1)=2*h2;a=-1*ones(1,n); b=4*ones(1,n);c=-1*ones(1,n); u=zeros(n+2,n+2); e=0.000000001;tic;for k=1:1000 er=0;t=u; for j=2:n+1 Ub=u(:,j); d(1:n)=f(2:n+1,j)+t(2:n+1,j-1)+u(2:n+1,j+1) ; % x=zg(a,b,c,d); % 用追趕法求解u(2:n+1,j)=x'

13、; er=er+abs(Ub-u(:,j); end if er/n2<e ,break;end % 判斷是否達(dá)到計(jì)算精度，如果達(dá)到則退出循環(huán)endT=tocksurf(u)function x=zg(a,b,c,d)n=length(b); %LU分解。u(1)=b(1); for k=2:n if u(k-1)=0 ,D=0,return; end l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*c(k-1);end% 追趕法求解之追過(guò)程，求解Ly=d。y(1)=d(1); for k=2:n y(k)=d(k)-l(k)*y(k-1);end if u(n)=0

14、,D=0,return; endx(n)=y(n)/u(n); for k=n-1:-1:1 x(k)=(y(k)-c(k)*x(k+1)/u(k); end結(jié)果：tyjak(8)T = 0.0703k = 120ans = 0 0 0 0 0 0 0 0 0 0 0 0.0300 0.0476 0.0577 0.0623 0.0623 0.0577 0.0476 0.0300 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0476 0 0 0.0577 0.0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0577 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0623 0.1045 0.1305 0.1429 0.1429 0.1305 0.1045 0.0623 0 0 0.0577 0.0961 0.1195 0.1305 0.1305 0.1195 0.0961 0.0577 0 0 0.0476 0.0780 0.0961 0.1045 0.1045 0.0961 0.0780 0.0

溫馨提示

  • 1. 本站所有資源如無(wú)特殊說(shuō)明，都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請(qǐng)下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等，如果需要附件，請(qǐng)聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙，網(wǎng)頁(yè)內(nèi)容里面會(huì)有圖紙預(yù)覽，若沒(méi)有圖紙預(yù)覽就沒(méi)有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫(kù)網(wǎng)僅提供信息存儲(chǔ)空間，僅對(duì)用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理，對(duì)用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯，并不能對(duì)任何下載內(nèi)容負(fù)責(zé)。
  • 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容，請(qǐng)與我們聯(lián)系，我們立即糾正。
  • 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時(shí)也不承擔(dān)用戶因使用這些下載資源對(duì)自己和他人造成任何形式的傷害或損失。

最新文檔

評(píng)論

0/150

提交評(píng)論