参考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 all;
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 = 10;%生成网格上的坐标
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 onfigure('name', 'abs err')
mesh(x,y, abs(u-u_analytical))
hold on
colorbar;
title('Absolute error')
hold on
