差分格式脚本文件:
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;
效果图: