Solving heat equation with Dirichlet boundary

Question

LM 2020-6-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/553285-solving-heat-equation-with-dirichlet-boundary

评论： darova 2020-6-24

在 MATLAB Online 中打开

Hey, I'm solving the heat equation

on a grid

for time

with inhomogeneous Dirichlet boundary conditions

. I'm using the implicit scheme for FDM, so I'm solving the Laplacian with the five-point-stencil, i.e.

where

are indices of the mesh. With the implicit scheme for the heat equation we get to solve

where A is the matrix representing the discretized Laplacian,

and F is zero if

is in the middle of the grid and

on the boundary because of the Dirichlet boundary which needs to be included in the scheme. My problem is, that the boundary of my plot is still equal to 0 and not

and I don't know why since I've changed the right hand side F according to the boundary. Even if I initialize U with

instead of

I get the boundary 0. Does somebody know why this happens?

Here is my code:

function diffusion
%heat equation is solved with lin. system of equations
T = 0.01;
J = 200; K = 200; L = (J-1)^2;
Delta_x = 1/J; Delta_t = 1/K;   %mesh-size, timestep-size
lambda = Delta_t/Delta_x^2;     %
end_time = floor(T/Delta_t);
R0 = 0.3;                       %radius of the circle
U = zeros(K+1,L);               
e = ones(J-1,1); E = ones(L,1);
X = spdiags([-e,4*e,-e],[-1,0,1],J-1,J-1);  %matrix to the left and right point in 5-point-stencil
A = sparse(L,L);                            %matrix for whole 5-point-stencil, it is a matrix which contains matrices
F = zeros(1,L);                             %right hand side
for j=1:J-2
    F((J-1)*j) = -1;                        %points on boundary of the grid, but not corner
    F((J-1)*j+1) = -1;
end
F(1) = -2; F(J-1) = -2; F((J-2)*(J-1)+1) = -2; F(L) = -2;   %points on each corner of the grid
for j = 1:(J-1):L
    A(j:j+J-2,j:j+J-2) = X;                 %puts X on the diagonal of A
end
A = A+spdiags([-E,-E],[-J+1,J-1],L,L);      %left and right we get -E because of upper and downer points of 5-point-stencil
B = speye(L)+lambda*A;
C = speye(L);
%initial condition
for j = 1:J-1
    for m =1:J-1
        U(1,j+(m-1)*(J-1)) = u_0([j,m]*Delta_x,R0);
    end
end
for k = 0:end_time
    show(U(k+1,:),J,Delta_x); drawnow;
    U(k+2,:) = B\(C*U(k+1,:)'+lambda*F');       %solve the implicit scheme
end
function val = u_0(x,R)
%initial condition is a circle
if ((x(1)-1/2)^2+(x(2)-1/2)^2)<=R^2
    u = 1;
else
    u = -1;
end
val = u;
function show(U,J,Delta_x)
U_mat = zeros(J+1,J+1); U_mat(2:J,2:J) = reshape(U,J-1,J-1)';
mesh(Delta_x*(0:J),Delta_x*(0:J),U_mat); %axis([0,1,0,1,0,1,0,1]); 