二維拋物線方程數(shù)值解法(ADI隱式交替法)方法

1、精選優(yōu)質(zhì)文檔-傾情為你奉上 ADI隱式交替法三種解法及誤差分析（一般的教材上只說第一種）理論部分參看孫志忠：偏微分方程數(shù)值解法注意：1. 最好不要直接看程序，中間很多公式很煩人的（一定要小心），我寫了兩天，終于寫對了。2. 中間：例如r*(u(i-1,m1,k)+u(i+1,m1,k)形式寫成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面會(huì)出錯(cuò)，我也不是很清楚為什么，可能由于舍入誤差，或者大數(shù)吃掉小數(shù)的影響。3. 下面有三個(gè)程序4. 具體理論看書，先仔細(xì)看書（孫志忠：偏微分方程數(shù)值解法）或者網(wǎng)上搜一些理論。Matlab程序：1.function u u0 p e x y t

2、=ADI1(h1,h2,m1,m2,n)%ADI解二維拋物線型偏微分方程(P-R交替隱式，截?cái)?%此程序用的是追趕法解線性方程組%h1為空間步長，h2為時(shí)間步長%m1,m2分別為x方向，y方向網(wǎng)格數(shù)，n為時(shí)間網(wǎng)格數(shù)%p為精確解，u為數(shù)值解,e為誤差%定義u0(i,j,k)=u(i,j,k+1/2),因?yàn)榫仃囍?，i,j,k必須全為整數(shù)x=(0:m1)*h1+0;%定義x0,y0,t0是為了f(x,t)=0的情況%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1 for i=1:m2+1 for j=1:m1+1 f(i,j,k)

3、=-1.5*exp(0.5*(x(j)+y(i)-t0(k); end endendfor i=1:m2+1 for j=1:m1+1 u(i,j,1)=exp(0.5*(x(j)+y(i); endendfor k=1:n+1 for i=1:m2+1 u(i,1 m1+1,k)=exp(0.5*y(i)-t(k) exp(0.5*(1+y(i)-t(k); u0(i,1 m1+1,k)=exp(0.5*y(i)-t0(k) exp(0.5*(1+y(i)-t0(k) ; endendfor k=1:n+1 for j=1:m1+1 u(1 m2+1,j,k)=exp(0.5*x(j)-t(

4、k) exp(0.5*(1+x(j)-t(k); u0(1 m2+1,j,k)=exp(0.5*x(j)-t0(k) exp(0.5*(1+x(j)-t0(k); endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);for k=1:n %外循環(huán)，先固定每一時(shí)間層，每一時(shí)間層上解一線性方程組% for i=2:m2 a=-r*ones(1,m1-1); c=a;a(1)=0;c(m1-1)=0; b=r2*ones(1,m1-1); d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k)+r1*u(i,2,k)+. h2*f(i,2,k)

5、; for l=2:m1-2 d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k)+r1*u(i,l+1,k)+. h2*f(i,l+1,k); %輸入部分系數(shù)矩陣，為0的矩陣元素不輸入%一定要注意輸入元素的正確性 end d(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k). +r1*u(i,m1,k)+h2*f(i,m1,k); for l=1:m1-2 %開始解線性方程組 消元過程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l);

6、 end u0(i,m1,k)=d(m1-1)/b(m1-1); %回代過程% for l=m1-2:-1:1 u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k)/b(l); end end for j=2:m1 a=-r*ones(1,m2-1); c=a;a(1)=0;c(m2-1)=0; b=r2*ones(1,m2-1); d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k)+r1*u0(2,j,k)+. h2*f(2,j,k); for l=2:m2-2 d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k)+r1

7、*u0(l+1,j,k)+. h2*f(l+1,j,k); %輸入部分系數(shù)矩陣，為0的矩陣元素不輸入%一定要注意輸入元素的正確性 end d(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k). +r1*u0(m2,j,k)+h2*f(m2,j,k); for l=1:m2-2 %開始解線性方程組 消元過程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l); end u(m2,j,k+1)=d(m2-1)/b(m2-1); %回代過程% for l

8、=m2-2:-1:1 u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1)/b(l); end endendfor k=1:n+1 for i=1:m2+1 for j=1:m1+1 p(i,j,k)=exp(0.5*(x(j)+y(i)-t(k); %p為精確解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k); %e為誤差 end endend2.function u p e x y t=ADI2(h1,h2,m1,m2,n)%ADI解二維拋物線型偏微分方程(D'Yakonov交替方向隱格式)%此程序用的是追趕法解線性方程組%h1為空間步長，h2為時(shí)

9、間步長%m1,m2分別為x方向，y方向網(wǎng)格數(shù)，n為時(shí)間網(wǎng)格數(shù)%p為精確解，u為數(shù)值解,e為誤差%定義u0(i,j,k)=u'(i,j,k)(引入的過渡層),因?yàn)榫仃囍?，i,j,k必須全為整數(shù)x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定義t0是為了f(x,y,t)=0的情況%for k=1:n+1 for i=1:m2+1 for j=1:m1+1 f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i)-t0(k); %編程時(shí)-t0(k)寫成了+t0(k),導(dǎo)致錯(cuò)誤; end endend%初始條件

10、for i=1:m2+1 for j=1:m1+1 u(i,j,1)=exp(0.5*(x(j)+y(i); endend%邊界條件for k=1:n+1 for i=1:m2+1 u(i,1 m1+1,k)=exp(0.5*y(i)-t(k) exp(0.5*(1+y(i)-t(k); endendr=h2/(h1*h1);r4=1+r;r5=r/2;for k=1:n for i=2:m2 u0(i,1 m1+1,k)=r4*u(i,1 m1+1,k+1)-r5*(u(i-1,1 m1+1,. k+1)+u(i+1,1 m1+1,k+1); endendfor k=1:n+1 for j=

11、1:m1+1 u(1 m2+1,j,k)=exp(0.5*x(j)-t(k) exp(0.5*(1+x(j)-t(k); endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;for k=1:n %外循環(huán)，先固定每一時(shí)間層，每一時(shí)間層上解一線性方程組% for i=2:m2 a=-r*ones(1,m1-1); c=a;a(1)=0;c(m1-1)=0; b=2*r4*ones(1,m1-1); d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+. u(i,3,k)+r2*u(i,2,k)+r3*(u(i-1,1

12、,k)+. u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k)+2*h2*f(i,2,k); for l=2:m1-2 d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+. u(i,l+2,k)+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+. u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k)+2*h2*f(i,l+1,k); %輸入部分系數(shù)矩陣，為0的矩陣元素不輸入%一定要注意輸入元素的正確性 end d(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+.

13、 u(i+1,m1,k)+u(i,m1+1,k)+r2*u(i,m1,k)+. r3*(u(i-1,m1-1,k)+. u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k)+2*h2*f(i,m1,k); for l=1:m1-2 %開始解線性方程組 消元過程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l); end %回代過程% u0(i,m1,k)=d(m1-1)/b(m1-1); for l=m1-2:-1:1 u0(i,l+1,k)=(d(l)-c(l)*u

14、0(i,l+2,k)/b(l); end end for j=2:m1 a=-r*ones(1,m2-1); c=a;a(1)=0;c(m2-1)=0; b=2*r4*ones(1,m2-1); d(1)=r*u(1,j,k+1)+2*u0(2,j,k); for l=2:m2-2 d(l)=2*u0(l+1,j,k); %輸入部分系數(shù)矩陣，為0的矩陣元素不輸入%一定要注意輸入元素的正確性 end d(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1); for l=1:m2-2 %開始解線性方程組 消元過程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1

15、)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l); end u(m2,j,k+1)=d(m2-1)/b(m2-1); %回代過程% for l=m2-2:-1:1 u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1)/b(l); end endendfor k=1:n+1 for i=1:m2+1 for j=1:m1+1 p(i,j,k)=exp(0.5*(x(j)+y(i)-t(k); %p為精確解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k); %e為誤差 end endend3.function u u0 p e x

16、y t=ADI5(h1,h2,m1,m2,n)%ADI解二維拋物線型偏微分方程(P-R交替隱式，未截?cái)?%此程序用的是追趕法解線性方程組%h1為空間步長，h2為時(shí)間步長%m1,m2分別為x方向，y方向網(wǎng)格數(shù)，n為時(shí)間網(wǎng)格數(shù)%p為精確解，u為數(shù)值解,e為誤差%定義u0(i,j,k)=u(i,j,k+1/2),因?yàn)榫仃囍?，i,j,k必須全為整數(shù)x=(0:m1)*h1+0;%定義x0,y0,t0是為了f(x,t)=0的情況%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1 for i=1:m2+1 for j=1:m1+1 f(i,

17、j,k)=-1.5*exp(0.5*(x(j)+y(i)-t0(k); end endendfor i=1:m2+1 for j=1:m1+1 u(i,j,1)=exp(0.5*(x(j)+y(i); endendfor k=1:n+1 for i=1:m2+1 u(i,1 m1+1,k)=exp(0.5*y(i)-t(k) exp(0.5*(1+y(i)-t(k); u1(i,1 m1+1,k)=exp(0.5*y(i)-t0(k) exp(0.5*(1+y(i)-t0(k) ; endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);for k=1:

18、n for i=2:m2 u0(i,1 m1+1,k)=u1(i,1 m1+1,k)-r2*(u(i-1,1 m1+1,k+1)-. 2*u(i,1 m1+1,k+1)+u(i+1,1 m1+1,k+1)-u(i-1,1 m1+1,k)+. 2*u(i,1 m1+1,k)-u(i+1,1 m1+1,k); endendfor k=1:n+1 for j=1:m1+1 u(1 m2+1,j,k)=exp(0.5*x(j)-t(k) exp(0.5*(1+x(j)-t(k); endendfor k=1:n %外循環(huán)，先固定每一時(shí)間層，每一時(shí)間層上解一線性方程組% for i=2:m2 a=-r*

19、ones(1,m1-1); c=a;a(1)=0;c(m1-1)=0; b=r3*ones(1,m1-1); d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k)+r1*u(i,2,k)+. h2*f(i,2,k); for l=2:m1-2 d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k)+r1*u(i,l+1,k)+. h2*f(i,l+1,k); %輸入部分系數(shù)矩陣，為0的矩陣元素不輸入%一定要注意輸入元素的正確性 end d(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k). +r1*u(i,

20、m1,k)+h2*f(i,m1,k); for l=1:m1-2 %開始解線性方程組 消元過程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l); end u0(i,m1,k)=d(m1-1)/b(m1-1); %回代過程% for l=m1-2:-1:1 u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k)/b(l); end end for j=2:m1 a=-r*ones(1,m2-1); c=a;a(1)=0;c(m2-1)=0; b=r3*ones(1,m2-1); d(

21、1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k)+r1*u0(2,j,k)+. h2*f(2,j,k); for l=2:m2-2 d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k)+r1*u0(l+1,j,k)+. h2*f(l+1,j,k); %輸入部分系數(shù)矩陣，為0的矩陣元素不輸入%一定要注意輸入元素的正確性 end d(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k). +r1*u0(m2,j,k)+h2*f(m2,j,k); for l=1:m2-2 %開始解線性方程組 消元過

22、程 a(l+1)=-a(l+1)/b(l); b(l+1)=b(l+1)+a(l+1)*c(l); d(l+1)=d(l+1)+a(l+1)*d(l); end u(m2,j,k+1)=d(m2-1)/b(m2-1); %回代過程% for l=m2-2:-1:1 u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1)/b(l); end endendfor k=1:n+1 for i=1:m2+1 for j=1:m1+1 p(i,j,k)=exp(0.5*(x(j)+y(i)-t(k); %p為精確解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k); %e為

23、誤差 end endend up e x y t=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001) t=1的誤差曲面下面是三種方法的誤差比較：1.u u0 p e x y t=ADI1(0.1,0.1,10,10,10)(P-R交替隱式，截?cái)?截?cái)嘀虚g過渡層用u(i,j,k+1/2）代替）（t=1時(shí)的誤差）2.u u0 p e x y t=ADI5(0.1,0.1,10,10,10)(P-R交替隱式，未截?cái)?(未截?cái)噙^渡層u(i,j,)=u(i,j,k+1/2)-h22/4*dy2dtu(i,j,k+1/2);)3.u p e x y t=A

24、DI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隱格式) surf(x,y,e(:,:,11)(表示t=1時(shí)的誤差）下面是相關(guān)數(shù)據(jù)：1: u u0 p e x y t=ADI1(0.1,0.1,10,10,10)(P-R交替隱式，截?cái)?截?cái)?中間過渡層用u(i,j,k+1/2)代替）e(:,:,11) = Columns 1 through 6 0 0 0 0 0 0 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0.

25、0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0 0 0 0 02.u u0 p e x y t=ADI5(0.1,0.1,10,10,10)(P-R交替隱式，未截?cái)?(未截?cái)噙^渡層u(i,j,)=u(i,j,k+1/2)-h22/4*dy2dtu(i,j,k+1/2);)e(:,:,11) = Columns 1 through 6 0 0 0 0 0 0 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0

26、 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0. 0. 0. 0. 0. 0 0 0 0 0 03.u p e x y t=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隱格式)e(:,:,11) = Columns 1 through 6 0 0 0 0 0 0 0 8.6469e-006 1.4412e-005 1.8364e-005 2.091e-005 2.2174e-005 0 1.4412e-005 2.4777e-005 3.2047e-005 3.6716e-005 3.8961e-0

27、05 0 1.8364e-005 3.2047e-005 4.1789e-005 4.8054e-005 5.1008e-005 0 2.091e-005 3.6716e-005 4.8054e-005 5.5353e-005 5.8764e-005 0 2.2174e-005 3.8961e-005 5.1008e-005 5.8764e-005 6.2389e-005 0 2.2118e-005 3.8698e-005 5.0523e-005 5.8126e-005 6.171e-005 0 2.055e-005 3.5581e-005 4.6157e-005 5.2942e-005 5.6197e-005 0 1.707e-005 2.8951e-005 3.7128e-005 4.2365e-005 4.4952e-005 0 1.0851e-005 1.7698e-005 2.2265e-005 2.5203e-005 2.672e-005 0 0 0 0 0 01.u u0 p e x y t=ADI1(0.1,0.1,10,10,10)(P-R交替隱式，截?cái)?截?cái)?中間過渡層用u(i,j,k+1/2)代替） Columns 7 through 11 0 0 0 0 0 0. 0. 0. 0. 0 0. 0. 0. 0.