matlab编写微分方程椭圆型方程（五点差分法）

$$
x=ih_1,

y=jh_2.
$$

% function chap2_fdm_elliptic_2D_nonzeroBD

%Poisson equation in 2D
%   -Delta u =—2*exp(x+y)，\Omega = (a, b)*(c, d)
%   u(x,0) = exp(x)，
%   u(x, 1) = exp(x+1)，
%   u(0,y) = exp(y)，
%   u(1, y) = exp(1+y)



clc,clear





global mLap M_V hx hy

a = 0;% lower bound of the domain
b = 1;% upper bound of the domain
c = 0;% lower bound of the domain
d = 1;% upper bound of the domain
M= 40;% number of subintervals of the partition of space
N = 30;
hx = (b-a)/M;
hy = (d-c)/N;
x = (a:hx: b)';
y = (c:hy: d)';
[X Y]= meshgrid(x, y);
x_in = (a+hx: hx: b-hx)';
y_in = (c+hy: hy: d-hy)';
[X_in Y_in] = meshgrid(x_in, y_in);

% the discretization of d ^2/dx^2
d_0 = -2*ones(M-1,1);
d_1 = 1*ones(M-1,1);
d_m1 = d_1;
d_xx = spdiags([d_m1 d_0 d_1],[-1 0 1],M-1,M-1);
d_xx = d_xx/hx^2;
I_Nm1 = speye(N-1,N-1);
D_xx = kron(I_Nm1, d_xx);

% the discretization of d 2/dy^2
d_0 = -2*ones(N-1,1);
d_1 = 1*ones(N-1,1);
d_m1 = d_1;
d_yy = spdiags([d_m1 d_0 d_1],[-1 0 1],N-1,N-1);
d_yy = d_yy/hy^2;
I_Mm1 = speye(M-1,M-1);
D_yy = kron(d_yy,I_Mm1);

% the discretization of Laplacian
mLap = -( D_xx + D_yy ) ;%29*39

% the right hand side with modification
F = f(X_in,Y_in);
F= F';
g_d = gx0(x_in);
g_u = gx1(x_in);
g_l = g0y(y_in);
g_r = g1y(y_in);
F(:,1)= F(:,1) + g_d/hy ^2;
F(:,N-1) = F(:,N-1) + g_u/hy^2;
F(1,:) = F(1,:) + g_l'/hx^2;
F(M-1,:)= F(M-1,:) + g_r'/hx^2;

% direct solver
%b=F^{~}
b = F(:);
%按列拉长，变成一个列向量


u = mLap\b;
% 矩阵除法求得u
U = reshape(u,M-1,N-1);%重排u顺序
U = U';

% plot the solution
u_e = u_exact(X,Y);
u_num = u_e;

%取精确值的边值，估算的内点u为U
u_num(2:N,2:M)= U;
figure(1)
mesh(X,Y,u_e);
figure(2)
mesh(X,Y,u_num);

% end

%——------——-——-——-—---------—-subroutines-----—---------
function y = f(x,y)
y =-2*exp(x+y);
end

function  u=gx0(x)
u = exp(x);
end

function  u=gx1(x)
u = exp(x+1);
end

function  u=g0y(y)
u = exp(y);
end

function  u=g1y(y)
u = exp(1+y);
end
%   u(x,0) = exp(x)，
%   u(x, 1) = exp(x+1)，
%   u(0,y) = exp(y)，
%   u(1, y) = exp(1+y)

%真解u = exp(x+y)
function u = u_exact(x,y)
u = exp(x+y);
end

注：第58行：

F= F';

转置的原因是M=40等分x，N=30等分y，$f(X_{i}n,Y_{i}n)$为$39\times 29$矩阵，需要转置一下将顺序放正。
举例来说：将x轴切为三份，y轴切为两份，得到矩阵应为$2\times 3$,但f（x，y）却为$3\times 2$，需将矩阵转置。

结果：

Figure1：

Figure2：