Matlab有限差分法求解狄利克雷(Dirichlet)边界的泊松(Poisson)问题，边界值为任意值

参考l链接：

有限差分法-二维泊松方程及其Matlab程序实现
弹性力学方程有限差分法matlab,泊松方程的有限差分法的MATLAB实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%          Matrix method for Poisson Equation         %%%%
%%% -[u_{xx} + u_{yy}]=f(x, y), xl < x < xr, yb < y < yt %%%%
%%%        u(x, y) = gl(xl,y) on left boundary,          %%%%  
%%%        u(x, y) = gr(xr,y) on left boundary,          %%%%  
%%%        u(x, y) = gb(x,yb) on left boundary,          %%%%  
%%%        u(x, y) = gt(x,yt) on left boundary,          %%%%  
%%%  Exact soln: u(x, y) = exp(pi*x)*sin(pi*y)           %%%%
%%%        Here f(x, y) = (pi^2-1)*exp(x)*sin(pi*y);     %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; 
clc;
close all;

fside = @(x, y) (pi^2-1)*exp(x).*sin(pi*y);
utrue = @(x, y) exp(x).*sin(pi*y);
uleft = @(x, y) exp(x).*sin(pi*y);
uright = @(x, y) exp(x).*sin(pi*y);
ubottom = @(x, y) exp(x).*sin(pi*y);
utop = @(x, y) exp(x).*sin(pi*y);

% 求解范围
xleft = 0.0;
xright = 2.0;
ybottom = 0.0;
ytop = 5;

%生成网格上的坐标
h=0.01; 
x=[xleft:h:xright]'; 
y=[ybottom:h:ytop]';
N=length(x)-1; 
M=length(y)-1;
[meshX,meshY]=meshgrid(x,y);

%解析解
u_analytical=exp(meshX).*sin(pi*meshY);

% 网格内部点(不包括边界)
meshX_in=meshX(2:M,2:N); 
meshY_in=meshY(2:M,2:N);

%在内部点上生成右端项f(x,y)
f=fside(meshX_in,meshY_in); 

% 左边界和右边界上的右端项
f(:,1)=f(:,1)+uleft(xleft,y(2:M))/h^2;
f(:,end)=f(:,end)+uright(xright,y(2:M))/h^2;

% 下边界和上边界上的右端项
temp2ub = ubottom(x(2:N), ybottom)/h^2;
temp2ut = utop(x(2:N),ytop)/h^2;
f(1,:)=f(1,:)+temp2ub';
f(end,:)=f(end,:)+temp2ut';

%构造矩阵D、C、A
I_element=ones(N-1,1);
C=1/h^2*spdiags([-I_element 4*I_element -I_element],[-1 0 1],N-1,N-1);
D=-1/h^2*eye(N-1);
I_mat=ones(M-1,1);
A=kron(eye(M-1),C)+kron(spdiags([I_mat I_mat],[-1 1],M-1,M-1),D);

%左除求解
f=f'; 
u=zeros(M+1,N+1);
u(2:M,2:N)=reshape(A\f(:),N-1,M-1)';  % 网格内部点上的解
u(:,1)=uleft(xleft,y);                % 左边界
u(:,end)=uleft(xright,y);             % 右边界

u(1,:)=ubottom(x,ybottom)';           % 左边界
u(end,:)=utop(x,ytop)';               % 右边界

%画图
figure('name','Exact')
mesh(x,y,u)
hold on
colorbar;
title('Exact solution')
hold on

figure('name', 'abs err')
mesh(x,y, abs(u-u_analytical))
hold on
colorbar;
title('Absolute error')
hold on