差分格式脚本文件：
1 tic; 2 clear 3 clc 4 M=32;%x的步数 5 N=16;%y的步数 6 h1=1/M;%x的步长 7 h2=1/N;%y的步长 8 x=0:h1:1; 9 y=0:h2:1; 10 u=zeros(M+1,N+1);%给数值解分配内存单元 11 U=u;%给精确解分配内存单元 12 u(1,:)=y.^3;%y边值 13 u(M+1,:)=1+y.^3;%y边值 14 u(:,1)=x.^3;%x边值。uo 15 u(:,N+1)=1+x.^3;%x边值。un 16 for i=1:M+1 17 for j=1:N+1 18 Accurate(i,j)=x(i)^3+y(j)^3;%精确解 19 end 20 end 21 fun=inline('-6*(x+y)','x','y'); 22 for i=1:M-1 23 for j=1:N-1 24 f(i,j)=fun(x(i+1),y(j+1)); 25 end 26 end 27 Numerical=u; 28 error=eye(M+1,N+1); 29 while norm(error,inf) >= 1e-5 30 for i=2:M 31 for j=2:N 32 Numerical(i,j)=(f(i-1,j-1)+N^2*u(i,j-1)+M^2*u(i-1,j)+M^2*u(i+1,j)+N^2*u(i,j+1))/(2*M^2+2*N^2); 33 end 34 end 35 error=Numerical-u; 36 u=Numerical; 37 end 38 [X,Y]=meshgrid(x,y); 39 subplot(1,3,1) 40 mesh(X,Y,Numerical'); 41 title('the image of Numerical solution') 42 xlabel('x');ylabel('y');zlabel('u'); 43 subplot(1,3,2) 44 mesh(X,Y,Accurate'); 45 title('the image of Accurate solution') 46 xlabel('x');ylabel('y');zlabel('U'); 47 subplot(1,3,3) 48 mesh(X,Y,(Numerical-Accurate)'); 49 title('the image of error solution') 50 xlabel('x');ylabel('y');zlabel('error'); 51 toc;
效果图：
紧差分格式：
1 tic; 2 clear 3 clc 4 M=100;%x的步数 5 N=100;%y的步数 6 h1=1/M;%x的步长 7 h2=1/N;%y的步长 8 x=0:h1:1; 9 y=0:h2:1; 10 u=zeros(M+1,N+1);%给数值解分配内存单元 11 U=u;%给精确解分配内存单元 12 u(1,:)=y.^3;%y边值 13 u(M+1,:)=1+y.^3;%y边值 14 u(:,1)=x.^3;%x边值。uo 15 u(:,N+1)=1+x.^3;%x边值。un 16 for i=1:M+1 17 for j=1:N+1 18 Accurate(i,j)=x(i)^3+y(j)^3;%精确解 19 end 20 end 21 b1=5/3*(M^2+N^2); 22 b2=1/6*(5*M^2-N^2); 23 b3=1/6*(5*N^2-M^2); 24 f=inline('-6*(x+y)','x','y'); 25 for i=1:M-1 26 for j=1:N-1 27 ABf(i,j)=(1/144)*(f(x(i),y(j))+10*f(x(i),y(j+1))+f(x(i),y(j+2))... 28 +10*f(x(i+1),y(j))+100*f(x(i+1),y(j+1))+10*f(x(i+1),y(j+2))... 29 +f(x(i+2),y(j))+10*f(x(i+2),y(j+1))+f(x(i+2),y(j+2))); 30 end 31 end 32 Numerical=u; 33 error=eye(M+1,N+1); 34 while norm(error,inf) >= 1e-10 35 for i=2:M 36 for j=2:N 37 Numerical(i,j)=(ABf(i-1,j-1)+1/20*b1*Numerical(i-1,j-1)+b3*Numerical(i,j-1)+1/20*b1*Numerical(i+1,j-1)... 38 +b2*Numerical(i-1,j)+b2*u(i+1,j)... 39 +1/20*b1*u(i-1,j+1)+b3*u(i,j+1)+1/20*b1*u(i+1,j+1))/b1; 40 end 41 end 42 error=Numerical-u; 43 u=Numerical; 44 end 45 [X,Y]=meshgrid(x,y); 46 subplot(1,3,1) 47 mesh(X,Y,Numerical'); 48 title('the image of Numerical solution') 49 xlabel('x');ylabel('y');zlabel('u'); 50 subplot(1,3,2) 51 mesh(X,Y,Accurate'); 52 title('the image of Accurate solution') 53 xlabel('x');ylabel('y');zlabel('U'); 54 subplot(1,3,3) 55 mesh(X,Y,(Numerical-Accurate)'); 56 title('the image of error solution') 57 xlabel('x');ylabel('y');zlabel('error'); 58 toc;
效果图：
