Laplace equation using central difference

12 次查看（过去 30 天）
Faraz Vossoughian
评论： darova 2020-5-4
Hi im trying to program a code that solves laplace equation using central difference. I already set up the boundary condition that corresponds to b matrix in the code, however im having hard time coding the A matrix that contains all the coefficients that looks like following. The A and b are found in the for loop in the code.
close all
h = 1.25; %Try different values & find the order of accuracy
% Geometry & boundary conditions
xmin = 0.0; xmax = 40.0;
ymin = 0.0; ymax = 40.0;
% Boundary conditions (In this problem, all boundaries are isothermal)
Tl = 75; %left
Tr = 50; %right
Tb = 0; %bottom
Tt = 100; %top
% Mesh (h is given as input)
nx = (xmax-xmin)/h - 1; %number of interior nodes in x-direction
ny = (ymax-ymin)/h - 1; %number of interior nodes in y-direction
N = nx*ny; %total number of interior nodes
% Allocate memory for A and b (sparse matrix & vector).
A = sparse(N,N); b = sparse(N,1);
% A=zeros(N,N);
% b=zeros(N,1);
% Loop through all the interior nodes, and fill A & b.
for k=1:N
A(k,k) = 4; % the diagonal element, corresponding to Tij, is always 4
%In the following, we look at the four neighbours of (i,j) that are
%involved in the 5-point formula
i = mod(k-1,nx)+1; j = (k-i)/nx+1; %get (i,j) from k.
disp(i)
disp(j)
if i==1
b(k,1)=Tl;
elseif i==nx
b(k,1)=Tr;
elseif j==1
b(k,1)=Tb;
elseif j==ny
b(k,1)=Tt;
end
end
% Solve Ax = b for nodal temperature values
T = full(A\b);
% Output Temperature at a sensor location (10 cm, 10 cm)
i = (0.25*(xmax-xmin))/h;
j = (0.25*(ymax-ymin))/h;
Tsensor = T(nx*(j-1)+i);
fprintf('h = %e, T(%d,%d) = %e.\n', h, i, j, Tsensor);
% Visualization
figure
[X, Y] = meshgrid(xmin:h:xmax, ymin:h:ymax);
Tvis = zeros(size(X)); %just for visualization
Tvis(:,1) = Tl; Tvis(:,nx+2) = Tr; Tvis(1,:) = Tb; Tvis(ny+2,:) = Tt;
for j = 2:nx+1
for i = 2:ny+1
Tvis(i,j) = T((i-2)*nx + j-1);
end
end
h = surf(X,Y,Tvis);
colormap('hot');
view(0,90);
shading interp
%set(h,'edgecolor','k');
cb = colorbar; cb.Label.String = 'Temperature (C)';
xlabel('X (cm)')
ylabel('Y (cm)')
zlabel('Temperature (C)')
daspect([1 1 1]);
  3 个评论
显示 1更早的评论
Faraz Vossoughian
Faraz Vossoughian 2020-5-4
Hi darobva, my question is set up differently, im really having a hard time codeing the coefficient matrix (A). could you please help
for k=1:N
A(k,k) = 4; % the diagonal element, corresponding to Tij, is always 4
% In the following, we look at the four neighbours of (i,j) that are
% involved in the 5-point formula
i = mod(k-1,nx)+1; j = (k-i)/nx+1; %get (i,j) from k.
disp('ass')
if i==1
b(k,1)=Tl;
elseif i==nx
b(k,1)=Tr;
elseif j==1
b(k,1)=Tb;
elseif j==ny
b(k,1)=Tt;
end
% if j<ny-2
% r=3+(j-1)
% A(i,nx-r)=-1
% end
%
if i>1 && j>1
if i==j
A(j-1,i)=-1;
A(j,i+1)=-1;
A(j,i-1)=-1;
A(j+1,i)=-1;
end
end
end
darova
darova 2020-5-4
Please look into the script i attached. Especially at these lines
for i = 2:m-1
for j = 2:n-1
ii = i + (j-1)*m;
A(ii,ii-1) = ci; % U(i-1,j)
A(ii,ii+1) = ci; % U(i+1,j)
A(ii,ii) = cij; % U(i,j)
A(ii,ii-m) = cj; % U(i,j-1)
A(ii,ii+m) = cj; % U(i,j+1)
end
end

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Geographic Plots 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by