二维泊松方程五点差分格式和九点差分格式的matlab求解

%二维泊松方程五点和九点差分格式
%-∇² = f(x,y)
%f(x,y) = pi^2/2*sin(pi*((x+1)/2))*sin(pi*((y+2)/2))
%真解 sin(pi*((x+1)/2)).*sin(pi*((y+2)/2))
%五点差分格式
%-Δₕu(i,j) = -[2*u(x,y)/h1^2+2*u(x,y)/h2^2-u(x,y-1)/h1^2-u(x,y+1)/h1^2-u(x-1,y)/h1^2-u(x+1,y)/h1^2] = fij
%九点差分格式
%-Δₕu(i,j)-1/12*(4*u(i,j)-2u(i-1,j)-2u(i+1,j)-2u(i,j+1)-2u(i,j-1)+u(i+1,j+1)+u(i+1,j-1)+u(i-1,j+1)+u(i-1,j-1))*(h1^2+h2^2)/(h1*h1*h2*h2) = fij+1/12*(h1^2*fxx(i,j)+h2^2*fyy(i,j))
%x,y的区间
a = -1;b = 1;%x
c = -2;d = 2;%y
%划分网格
h = 0.05;%定步长
N1 = (b-a)/h;
N2 = (d-c)/h;
%N1 = 300;%指定的网格x轴被划分多少份
%N2 = 200;%指定的网格y轴被划分多少份
h1 = (b-a)/N1;
h2 = (d-c)/N2;
x = a+[1:N1-1]*h1;
y = c+[1:N2-1]*h2;
%五点差分矩阵
%五点差分格式系数放到对角线上
e = ones((N1-1)*(N2-1),1);
A = spdiags([e/(h1^2),e/(h2^2),-2*e/(h1^2)-2*e/(h2^2),e/(h2^2),e/(h1^2)], [-(N2-1),-1,0,1,N2-1] ,(N1-1)*(N2-1),(N1-1)*(N2-1));
%处理边界
% 假设N1=N2=5,A
% 2  1  0  0  1  0  0  0  0  0  0  0  0  0  0  0
% 1  2  1  0  0  1  0  0  0  0  0  0  0  0  0  0
% 0  1  2  1  0  0  1  0  0  0  0  0  0  0  0  0
% 0  0  1  2 #1  0  0  1  0  0  0  0  0  0  0  0

% 1  0  0 ##1 2  1  0  0  1  0  0  0  0  0  0  0
% 0  1  0  0  1  2  1  0  0  1  0  0  0  0  0  0
% 0  0  1  0  0  1  2  1  0  0  1  0  0  0  0  0
% 0  0  0  1  0  0  1  2##1  0  0  1  0  0  0  0

% 0  0  0  0  1  0  0##1  2  1  0  0  1  0  0  0
% 0  0  0  0  0  1  0  0  1  2  1  0  0  1  0  0
% 0  0  0  0  0  0  1  0  0  1  2  1  0  0  1  0
% 0  0  0  0  0  0  0  1  0  0  1  2##1  0  0  1

% 0  0  0  0  0  0  0  0  1  0  0  1  2  1  0  0
% 0  0  0  0  0  0  0  0  0  1  0  0  1  2  1  0
% 0  0  0  0  0  0  0  0  0  0  1  0  0  1  2  1
% 0  0  0  0  0  0  0  0  0  0  0 #1  0  0  1  2

%将#1处替换为0
A(N2-1,N2) = 0;
A((N1-2)*(N2-1)+1,(N1-2)*(N2-1))=0;
%将##1处替换为0
for i = 1:N1-1-2
    A((N2-1)*i+1,(N2-1)*i) = 0;
    A((N2-1)*(i+1),(N2-1)*(i+1)+1) = 0;
end

%等式右端f(x,y)赋值
f = @(x, y) sin(pi*((x+1)/2))*sin(pi*((y+2)/2))*pi*pi/2;
F = zeros((N1-1)*(N2-1),1);
n=1;
for i = 1:N1-1
    for j = 1:N2-1
        F(n) = f(a+i*h1,c+j*h2);
        n = n+1;
    end
end

%五点差分格式求解
A = -1*A;
u=A\F;
%五点格式绘图
[X,Y]=meshgrid(x,y);
U = reshape(u,N2-1,N1-1);
%求解图
figure(1)
surf(X,Y,U)
xlabel('x');
ylabel('y');
zlabel('五点格式U(x,y)');
axis([a b c d -2 2])
%真解图
figure(2)
Z = sin(pi*((X+1)/2)).*sin(pi*((Y+2)/2));
surf(X, Y, Z); % 绘制三维图像
xlabel('x');
ylabel('y');
zlabel('真解');
axis([a b c d -2 2])
%error图
figure(3)
P = Z-U;
surf(X, Y, P);
xlabel('x');
ylabel('y');
zlabel('五点格式error');

%九点差分格式矩阵
e1 = ones(N1-1,1);
e2 = ones(N2-1,1);
C = spdiags([e1,-2*e1,e1],[-1,0,1],N1-1,N1-1);
D = spdiags([e2,-2*e2,e2],[-1,0,1],N2-1,N2-1);
E = kron(C,D);
%假设N1=N2=5,E
%  4  -2   0   0  -2   1   0   0   0   0   0   0   0   0   0   0
% -2   4  -2   0   1  -2   1   0   0   0   0   0   0   0   0   0
%  0  -2   4  -2   0   1  -2   1   0   0   0   0   0   0   0   0
%  0   0  -2   4   0   0   1  -2   0   0   0   0   0   0   0   0

% -2   1   0   0   4  -2   0   0  -2   1   0   0   0   0   0   0
%  1  -2   1   0  -2   4  -2   0   1  -2   1   0   0   0   0   0
%  0   1  -2   1   0  -2   4  -2   0   1  -2   1   0   0   0   0
%  0   0   1  -2   0   0  -2   4   0   0   1  -2   0   0   0   0

%  0   0   0   0  -2   1   0   0   4  -2   0   0  -2   1   0   0
%  0   0   0   0   1  -2   1   0  -2   4  -2   0   1  -2   1   0
%  0   0   0   0   0   1  -2   1   0  -2   4  -2   0   1  -2   1
%  0   0   0   0   0   0   1  -2   0   0  -2   4   0   0   1  -2

%  0   0   0   0   0   0   0   0  -2   1   0   0   4  -2   0   0
%  0   0   0   0   0   0   0   0   1  -2   1   0  -2   4  -2   0
%  0   0   0   0   0   0   0   0   0   1  -2   1   0  -2   4  -2
%  0   0   0   0   0   0   0   0   0   0   1  -2   0   0  -2   4
% E = E+1-1;
%等式右端fxx+fyy赋值
syms x1 x2 
f0 = sin(pi*((x1+1)/2))*sin(pi*((x2+2)/2))*pi*pi/2;
f1 = diff(f0,x1,2);
f2 = diff(f0,x2,2);
f3 = matlabFunction(f1);
f4 = matlabFunction(f2);
n=1;
for i = 1:N1-1
    for j = 1:N2-1
        F(n) = F(n)+(h1*h1*f3(a+i*h1,c+j*h2)+h2*h2*f4(a+i*h1,c+j*h2))/12;
        n = n+1;
    end
end
%九点差分格式求解
A1 = -1*A-(h1^2+h2^2)/(h1*h1*h2*h2)*(1/12)*E;
u1=A\F;
%九点格式绘图
[X,Y]=meshgrid(x,y);
U1 = reshape(u1,N2-1,N1-1);
%求解图
figure(4)
surf(X,Y,U1)
xlabel('x');
ylabel('y');
zlabel('九点格式U(x,y)');
axis([a b c d -2 2])
%真解图
figure(5)
Z = sin(pi*((X+1)/2)).*sin(pi*((Y+2)/2));
surf(X, Y, Z); % 绘制三维图像
xlabel('x');
ylabel('y');
zlabel('真解');
axis([a b c d -2 2])
%error图
figure(6)
P1 = Z-U1;
surf(X, Y, P1);
xlabel('x');
ylabel('y');
zlabel('九点格式error');