Mtlab:抛物型方程的交替方向隐格式(ADI)
1 tic; 2 clear 3 clc 4 M=[5,10,20,40,80]; 5 N=M; 6 for p=1:length(M) 7 h=1/M(p);% 这里定义空间步长等距 8 tau=1/N(p); % 时间步长 9 x=0:h:1; 10 y=0:h:1; 11 t=0:tau:1; 12 Numerical(M(p)+1,M(p)+1,N(p)+1)=0;%u 13 numerical(M(p)+1,M(p)-1)=0;%u* 14 %-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 15 % 求解u*ij和uij过程中构造三对角矩阵 16 % a 表示下对角线元素 17 % b 表示主对角线元素 18 % c 表示上对角线元素 19 a=-tau/(2*h^2)*ones(M(p)-2,1); 20 b=(tau/h^2+1)*ones(M(p)-1,1); 21 c=-tau/(2*h^2)*ones(M(p)-2,1); 22 %-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* 23 for i=1:M(p)+1 24 for j=1:M(p)+1 25 Numerical(i,j,1)=exp(1/2*(x(i)+y(j)));% 初值Numerical（x,y,0）=u(i,j,0) 26 end 27 end 28 for j=1:M(p)+1 29 for k=1:N(p)+1 30 Numerical(1,j,k)=exp(1/2*y(j)-t(k));% 边值Numerical（0,y,t）=u(0,j,k) 31 end 32 end 33 for j=1:M(p)+1 34 for k=1:N(p)+1 35 Numerical(M(p)+1,j,k)=exp(1/2*(1+y(j))-t(k));% 边值Numerical（1,y,t）=u(m,j,k) 36 end 37 end 38 for i=1:M(p)+1 39 for k=1:N(p)+1 40 Numerical(i,1,k)=exp(1/2*x(i)-t(k));% 边值Numerical（x,0,t）=u(i,0,k) 41 end 42 end 43 for i=1:M(p)+1 44 for k=1:N(p)+1 45 Numerical(i,M(p)+1,k)=exp(1/2*(1+x(i))-t(k));% 边值Numerical（x,1,t）=u(i,m,k) 46 end 47 end 48 f=inline('-3/2*exp(1/2*(x+y)-t)','x','y','t'); 49 fun=inline('exp(1/2*(x+y)-t)','x','y','t'); 50 for i=1:M(p)+1 51 for j=1:M(p)+1 52 for k=1:M(p)+1 53 Accurate(i,j,k)=fun(x(i),y(j),t(k)); 54 end 55 end 56 end 57 for k=1:N(p); 58 for j=1:M(p)-1;% 固定j 59 numerical(1,j)=-tau/(2*h^2)*Numerical(1,j,k+1)+(tau/h^2+1)*Numerical(1,j+1,k+1)-tau/(2*h^2)*Numerical(1,j+2,k+1);% u*0j 60 numerical(M(p)+1,j)=-tau/(2*h^2)*Numerical(M(p)+1,j,k+1)+(tau/h^2+1)*Numerical(M(p)+1,j+1,k+1)-tau/(2*h^2)*Numerical(M(p)+1,j+2,k+1);% u*mj 61 for i=1:M(p)-1 62 numerical_right_vector(i,1)=tau*f(x(i+1),y(j+1),t(k)+tau/2)+Numerical(i+1,j+1,k)... 63 +tau/(2*h^2)*(Numerical(i,j+1,k)-2*Numerical(i+1,j+1,k)+Numerical(i+2,j+1,k))... 64 +tau/(2*h^2)*(Numerical(i+1,j,k)-2*Numerical(i+1,j+1,k)+Numerical(i+1,j+2,k))... 65 +tau^2/(4*h^4)*(Numerical(i,j,k)-2*Numerical(i+1,j,k)+Numerical(i+2,j,k))... 66 +tau^2/(4*h^4)*(-2*Numerical(i,j+1,k)+4*Numerical(i+1,j+1,k)-2*Numerical(i+2,j+1,k))... 67 +tau^2/(4*h^4)*(Numerical(i,j+2,k)-2*Numerical(i+1,j+2,k)+Numerical(i+2,j+2,k)); 68 end 69 numerical_right_vector(1,1)=numerical_right_vector(1,1)+tau/(2*h^2)*numerical(1,j); 70 numerical_right_vector(M(p)-1,1)=numerical_right_vector(M(p)-1,1)+tau/(2*h^2)*numerical(M(p)+1,j); 71 numerical(2:M(p),j)=chase(a,b,c,numerical_right_vector); 72 end 73 for i=1:M(p)-1 % 固定i 74 for j=1:M(p)-1 75 Numerical_right_vector(j,1)=numerical(i+1,j); 76 end 77 Numerical_right_vector(1,1)=Numerical_right_vector(1,1)+tau/(2*h^2)*Numerical(i+1,1,k+1); 78 Numerical_right_vector(M(p)-1,1)=Numerical_right_vector(M(p)-1,1)+tau/(2*h^2)*Numerical(i+1,M(p)+1,k+1); 79 Numerical(i+1,2:M(p),k+1)=chase(a,b,c,Numerical_right_vector); 80 end 81 end 82 error=abs(Numerical(:,:,M(p)+1)-Accurate(:,:,M(p)+1)); 83 error_inf(p)=max(max(error)); 84 figure(p) 85 [X,Y]=meshgrid(y,x); 86 subplot(1,3,1) 87 surf(X,Y,Accurate(:,:,M(p))); 88 xlabel('x');ylabel('y');zlabel('Numerical'); 89 title('the image of Accurate rusult'); 90 grid on; 91 subplot(1,3,2) 92 surf(X,Y,Numerical(:,:,M(p))); 93 xlabel('x');ylabel('y');zlabel('Numerical'); 94 title('the image of Numerical'); 95 grid on; 96 subplot(1,3,3) 97 surf(X,Y,error); 98 xlabel('x');ylabel('y');zlabel('error'); 99 title('the image of error Numerical'); 100 grid on; 101 end 102 for k=2:length(M) 103 X=error_inf(k-1)/error_inf(k); 104 Norm(k-1)=log2(X); 105 end 106 figure(length(N)+1) 107 plot(1:length(N)-1,Norm,'-b^'); 108 xlabel('序号');ylabel('误差阶数'); 109 title('ADI格式误差阶'); 110 grid on; 111 toc;
取t=1，N=5,10,20,40,80；真解与数值解的结果图：
N=5:
N=10;
N=20;
N=40;
N=80;
误差阶数：
胡冬冬
