Solving Laplace equation using Finite difference
57 ビュー (過去 30 日間)
古いコメントを表示
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')
0 件のコメント
採用された回答
Torsten
2022 年 11 月 28 日
編集済み: Torsten
2022 年 11 月 28 日
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.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Boundary Conditions についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!