ADI.rar_ADI隐式交替法_P-R_交替方向法_交替隐式_隐式差分 matlab

clear;format long r=1;aa=1/(2*pi^2); n=19;M=((n+1)^2)/r; h=1/(n+1); x=0:h:1;y=0:h:1;t=0; r=r*aa; A=zeros(n,n);a=1+r;b=-r/2; F=zeros(n,1);err=1; U=zeros(n+2,n+2); U_exact=zeros(n+2,n+2); %输入稀疏矩阵A for i=1:n-1 A(i,i)=a; A(i+1,i)=b; A(i,i+1)=b; end A(n,n)=a; %消元法化简矩阵A for k=1:n-1 A(k+1,k)=A(k+1,k)/A(k,k); A(k+1,k+1)=A(k+1,k+1)-A(k,k+1)*A(k+1,k); end %initial value for i=1:n+2 for j=1:n+2 U(i,j)=sin(pi*x(i))*sin(pi*y(j)); U_exact(i,j)=exp(-1)*sin(pi*x(i))*sin(pi*y(j)); end end %ADI for k=1:M for j=2:n+1 %deal with the right const for i=2:n+1 F(i-1)=r/2*U(i,j+1)+(1-r)*U(i,j)+r/2*U(i,j-1); end for i=2:n F(i)=F(i)-F(i-1)*A(i,i-1); end %huidai U(n+1,j)=F(n)/A(n,n); for i=n:-1:2 U(i,j)=(F(i-1)-A(i-1,i)*U(i+1,j))/A(i-1,i-1); end end for i=2:n+1 %deal with the right const for j=2:n+1 F(j-1)=r/2*U(i-1,j)+(1-r)*U(i,j)+r/2*U(i+1,j); end for j=2:n F(j)=F(j)-F(j-1)*A(j,j-1); end %huidai U(i,n+1)=F(n)/A(n,n); for j=n:-1:2 U(i,j)=(F(j-1)-A(j-1,j)*U(i,j+1))/A(j-1,j-1); end end end %error err1=0; for i=1:n+2 %计算误差 for j=1:n+2 err1=err1+(U(i,j)-U_exact(i,j))^2; end end err2=sqrt(err1)*h figure surf(x,y,U)%画近似值函数图像 xlabel('x') ylabel('y') zlabel('U') title('n=9时近似解的函数图像') figure surfc(x,y,U_exact) %画精确值函数图像 xlabel('x') ylabel('y') zlabel('U_exact') title('n=9时精确解的函数图像')