Solving Laplace equation using Finite difference

50 views (last 30 days)
Greetings all,
I'm trying to solve the following problem using a finite differnce iterative scheme. I wrote the following code which seems to give me a solution that does not vary with changing the length. I end up getting the same solution regardless of the domain length which contradicts the solution from the pde tool. My code is below, your help is appericiated. I think the problem is in how I implement the loop maybe, I'm not sure.
clear all;
clc;
% Apply user input here
L = 3;
H = 1;
pts = 80;
% In case of undivisible by 3 grid
int_part = fix(pts/3);
rem_part = pts - (int_part*3);
% Grid and tolerance
tol = 1.0e-5;
if rem_part ~= 0
x1 = int_part;
x2 = (x1 * 2);
x3 = (x1 * 3) + rem_part;
else
x1 = int_part;
x2 = x1*2;
x3 = x1*3;
end
x = 0:L/(pts-1):L;
y = 0:H/(pts-1):H;
%[X,Y] = meshgrid(x,y);
%% Initialization
U = zeros(pts);
U0 = zeros(pts);
eta = 1; % dummy
%% boundary conditions
r1 = 2:x1;
r2 = (x1+1):x2;
r3 = (x2+1):(x3-1);
rx = 2:pts-1;
U(pts,1:x1) = -100;
U(1,r2) = 100;
U(pts,(x2+1):(x3)) = -100;
%% Solution loop
i=2:pts-1;
j=2:pts-1;
while eta > tol
U(i,j) = 0.25*(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)); % Inner Domain
U(rx,1) = 0.25*(U(rx,1)+U(rx,1)+U(rx,2)+U(rx,2)); % b.c -> ux = 0
U(rx,pts) = 0.25*(U(rx,pts)+U(rx,pts)+U(rx,pts-1)+U(rx,pts-1)); % b.c -> ux = 0
U(1,r1) = 0.25*(U(2,r1)+U(2,r1)+U(1,r1-1)+U(1,r1+1)); % b.c -> uy = 0
U(1,r3) = 0.25*(U(2,r3)+U(2,r3)+U(1,r3-1)+U(1,r3+1)); % b.c -> uy = 0
U(pts,r2) = 0.25*(U(pts-1,r2)+U(pts-1,r2)+U(pts,r2-1)+U(pts,r2+1)); % b.c -> uy = 0
U(1,1) = 0.25*(U(2,1)+U(2,1)+U(1,2)+U(1,2)); % corner point
U(1,pts) = 0.25*(U(2,pts)+U(2,pts)+U(1,pts-1)+U(1,pts-1)); % corner point
eta = max(max(abs(U-U0)));
U0 = U;
end
contourf(x,y,U)
colorbar
colormap('jet')

Accepted Answer

Torsten
Torsten on 28 Nov 2022
Edited: Torsten on 28 Nov 2022
If L = 3 and H = 1, you cannot choose the same number of points in x and y direction if you update U as if you use a grid with equal cell sizes in x and y direction.
So either you choose nx = 3*ny for nx, ny being the number of cells (not points) in x and y direction, respectively, or you update your U values according to
(U(i+1,j)-2*U(i,j)+U(i-1,j))/dx^2 + (U(i,j+1)-2*U(i,j)+U(i,j-1))/dy^2 = 0
thus
U(i,j) = 1/(2/dx^2+2/dy^2) * ( (U(i+1,j)+U(i-1,j))/dx^2 + (U(i,j+1)+U(i,j-1))/dy^2 )
with
dx = L/nx, dy = H/ny
(Similar for the points where derivatives are prescribed).
And of course all the U's on the right-hand side of your updates
U(i,j) = 0.25*(U(i-1,j)+U(i+1,j)+U(i,j-1)+U(i,j+1)); % Inner Domain
U(rx,1) = 0.25*(U(rx,1)+U(rx,1)+U(rx,2)+U(rx,2)); % b.c -> ux = 0
U(rx,pts) = 0.25*(U(rx,pts)+U(rx,pts)+U(rx,pts-1)+U(rx,pts-1)); % b.c -> ux = 0
U(1,r1) = 0.25*(U(2,r1)+U(2,r1)+U(1,r1-1)+U(1,r1+1)); % b.c -> uy = 0
U(1,r3) = 0.25*(U(2,r3)+U(2,r3)+U(1,r3-1)+U(1,r3+1)); % b.c -> uy = 0
U(pts,r2) = 0.25*(U(pts-1,r2)+U(pts-1,r2)+U(pts,r2-1)+U(pts,r2+1)); % b.c -> uy = 0
U(1,1) = 0.25*(U(2,1)+U(2,1)+U(1,2)+U(1,2)); % corner point
U(1,pts) = 0.25*(U(2,pts)+U(2,pts)+U(1,pts-1)+U(1,pts-1)); % corner point
must be changed to U0's.
  1 Comment
Khaled Mohamed
Khaled Mohamed on 28 Nov 2022
Thank you. That was the problem. I was confused. Your help is really apperciated.

Sign in to comment.

More Answers (0)

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!