# Trying to solve 2 dimensional Partial differential equation using Finite Difference Method

조회 수: 91 (최근 30일)
Tristofani Agasta . 2022년 3월 21일
편집: Torsten . 2023년 9월 28일
Currently I study about finite difference for 1d and 2d partial differential equation. I finish my code by trying to follow the algorithm my lecturer gave to me. The difference is, I add some conditional for some nodes which are located at boundaries (at the top and the right where the value supposedly be 1, not 0). But why my graph seems wrong? This method seems similar to Sandip Mazumder book and Youtube tutorial.
clc; clear all; close all
Ny=30;Nx=30;
dx=0.01;dy=0.01;
xa=0:dx:(Nx-1)*dx;
ya=0:dy:(Ny-1)*dy;
yb=(Ny-1)*dy:-dy:0;
xb=(Nx-1)*dx:-dx:0;
a=1/(dx)^2; c=1/(dy^2);
b=-2*(a+c);
%Create matrices A, B and solution
[A,B]= matriks(a, b,c, Nx, Ny);
solx =inv(A)*B;
for ii=1:Ny;
for jj=1:Nx;
k=(jj-1)*Ny+ii;
sol(ii, jj)=solx(k);
end
end
%Showing the graph
[X, Y] = meshgrid(xb, yb);
surface(X, Y, sol); colormap
shading interp; axis ('equal')
xlim([0 max(xa)]);ylim([0 max(ya)])
xlabel('Sumbu X'); ylabel('Sumbu Y')
function [A,B]=matriks(a,b,c,Nx,Ny)
B=zeros(Nx*Ny, 1);
A=eye(Nx*Ny);
dx=0.01
x=0:dx:1
y=0:dx:1
for ii=1:Nx;
for jj=1:Ny
if (ii>1) && (ii<Nx) && (jj>1) && (jj<Ny) % Insides
k=(jj-1)*Ny+ii;
B(k,1)=0;
A(k,k)=b;
A(k, k-1)=a;A(k,k+1)=a;
A(k, k-Ny)=c;A(k, k+Ny)=c;
elseif (jj==Ny) && (ii>1) && (ii<Nx) % Top boundary
k=(jj-1)*Ny+ii;
B(k,1)=y(jj);
A(k,k)=b;
A(k, k-1)=a;A(k,k+1)=a;
A(k, k-Ny)=c;
elseif ( ii==Nx ) &&( jj<Ny) && ( jj > 1 ) % Right boundary
k=( jj - 1 )*Ny+ii;
B( k , 1 )=x(jj);
A( k , k )=b;
A( k, k - 1 )=a;
if k < (Ny*Nx)-Ny
A( k , k - Ny )=c;A(k, k+Ny)=c;
end
end
end
end
end

댓글을 달려면 로그인하십시오.

### 답변 (2개)

Torsten 2022년 3월 21일
편집: Torsten 님. 2022년 3월 21일
For dx = dy = 0.01, Nx = Ny = 101, not 30 in your code. I just realized this after setting up the code below.
dx = 0.01;
dy = 0.01;
x = 0:dx:1;
y = 0:dy:1;
nx = numel(x);
ny = numel(y);
a =1/dx^2;
c =1/dy^2;
b =-2*(a+c);
A = zeros(nx*ny,nx*ny);
B = zeros(nx*ny,1);
% Boundaries
% Boundary values at y = 0
for ix = 1:nx
A(ix,ix) = 1.0;
B(ix) = 0.0;
end
% Boundary values at x = 0
for iy = 2:ny-1
k = nx*(iy-1) + 1;
A(k,k) = 1.0;
B(k) = 0.0;
end
% Boundary values at x = 1
for iy = 2:ny-1
k = nx*iy;
A(k,k) = 1.0;
B(k) = y(iy);
end
% Boundary values at y = 1
for ix = 1:nx
k = nx*(ny-1) + ix;
A(k,k) = 1.0;
B(k) = x(ix);
end
% Inner grid points
for iy = 2:ny-1
for ix = 2:nx-1
k = (iy-1)*nx + ix;
A(k,k) = b;
A(k,k+1) = a;
A(k,k-1) = a;
A(k,k+nx) = c;
A(k,k-nx) = c;
end
end
u = A\B;
%u
U = zeros(nx,ny);
for iy = 1:ny
for ix = 1:nx
k = (iy-1)*nx + ix;
U(ix,iy) = u(k);
end
end
[X,Y] = meshgrid(x,y);
surf(X, Y, U);
##### 댓글 수: 5이전 댓글 3개 표시이전 댓글 3개 숨기기
Torsten 2022년 3월 22일
편집: Torsten 님. 2022년 3월 22일
Which iterations do you mean ? The loops to set up A and B ?
Tristofani Agasta 2022년 3월 23일
Both A and B, I thought I can put all nodes on each boundary in same 'for' loop. I never thought I can use for ii and for jj for each boundary separately, with different k index

댓글을 달려면 로그인하십시오.

Torsten 2023년 9월 27일
For those who might be interested in a finite volume code for the equation above.
Skerdi Hymeraj asked for such code, but deleted the given answer.
dx = 0.01;
dy = 0.01;
x = dx/2:dx:1-dx/2;
y = dy/2:dy:1-dy/2;
nx = numel(x);
ny = numel(y);
%A = zeros(nx*ny);
b = zeros(nx*ny,1);
index = 0;
% Points in contact to boundaries
% i = 1, j = 1
k = 1;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = ( -dy/dx - dy/(dx/2) - dx/dy - dx/(dy/2) );
index = index + 1;
irc(index) = k;
icc(index) = k+1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k+nx;
mat(index) = dx/dy;
%A(k,k) = ( -dy/dx - dy/(dx/2) - dx/dy - dx/(dy/2) );
%A(k,k+1) = dy/dx;
%A(k,k+nx) = dx/dy;
b(k) = -dy/(dx/2) * bcfun(0,y(1)) -dx/(dy/2) * bcfun(x(1),0);
% i = nx, j = 1
k = nx;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = ( -dy/(dx/2) - dy/dx - dx/dy - dx/(dy/2) );
index = index + 1;
irc(index) = k;
icc(index) = k-1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k + nx;
mat(index) = dx/dy;
%A(k,k) = ( -dy/(dx/2) - dy/dx - dx/dy - dx/(dy/2) );
%A(k,k-1) = dy/dx;
%A(k,k+nx) = dx/dy;
b(k) = -dy/(dx/2) * bcfun(1,y(1)) -dx/(dy/2) * bcfun(x(nx),0);
% i = 1, j = ny
k = (ny-1)*nx+1;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = (-dy/dx - dy/(dx/2) - dx/(dy/2) - dx/dy);
index = index + 1;
irc(index) = k;
icc(index) = k+1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k-nx;
mat(index) = dx/dy;
%A(k,k) = (-dy/dx - dy/(dx/2) - dx/(dy/2) - dx/dy);
%A(k,k+1) = dy/dx;
%A(k,k-nx) = dx/dy;
b(k) = -dy/(dx/2) * bcfun(0,y(ny)) -dx/(dy/2) * bcfun(x(1),1);
% i = nx, j = ny
k = nx*ny;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = (-dy/(dx/2) - dy/dx - dx/(dy/2) - dx/dy);
index = index + 1;
irc(index) = k;
icc(index) = k-1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k-nx;
mat(index) = dx/dy;
%A(k,k) = (-dy/(dx/2) - dy/dx - dx/(dy/2) - dx/dy);
%A(k,k-1) = dy/dx;
%A(k,k-nx) = dx/dy;
b(k) = -dy/(dx/2) * bcfun(1,y(ny)) -dx/(dy/2) * bcfun(x(nx),1);
% 1 < i < nx, j = 1
for ix = 2:nx-1
k = ix;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = -dy/dx - dy/dx - dx/dy -dx/(dy/2);
index = index + 1;
irc(index) = k;
icc(index) = k-1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k+1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k+nx;
mat(index) = dx/dy;
%A(k,k) = -dy/dx - dy/dx - dx/dy -dx/(dy/2);
%A(k,k-1) = dy/dx;
%A(k,k+1) = dy/dx;
%A(k,k+nx) = dx/dy;
b(k) = -dx/(dy/2) * bcfun(x(ix),0);
end
% 1 < i < nx, j = ny
for ix = 2:nx-1
k = (ny-1)*nx + ix;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = -dy/dx -dy/dx -dx/(dy/2) -dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k-1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k+1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k-nx;
mat(index) = dx/dy;
%A(k,k) = -dy/dx -dy/dx -dx/(dy/2) -dx/dy;
%A(k,k-1) = dy/dx;
%A(k,k+1) = dy/dx;
%A(k,k-nx) = dx/dy;
b(k) = -dx/(dy/2) * bcfun(x(ix),1);
end
% i = 1, 1 < j < ny
for iy = 2:ny-1
k = 1 + (iy-1)*nx;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = -dy/dx - dy/(dx/2) - dx/dy - dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k+1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k-nx;
mat(index) = dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k+nx;
mat(index) = dx/dy;
%A(k,k) = -dy/dx - dy/(dx/2) - dx/dy - dx/dy;
%A(k,k+1) = dy/dx;
%A(k,k-nx) = dx/dy;
%A(k,k+nx) = dx/dy;
b(k) = -dy/(dx/2) * bcfun(0,y(iy));
end
% i = nx, 1 < j < ny
for iy = 2:ny-1
k = nx*iy;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = -dy/(dx/2) - dy/dx - dx/dy - dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k-1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k-nx;
mat(index) = dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k+nx;
mat(index) = dx/dy;
%A(k,k) = -dy/(dx/2) - dy/dx - dx/dy - dx/dy;
%A(k,k-1) = dy/dx;
%A(k,k-nx) = dx/dy;
%A(k,k+nx) = dx/dy;
b(k) = -dy/(dx/2) * bcfun(1,y(iy));
end
% Inner grid points
% 1 < ix < nx, 1 < iy < ny
for ix = 2:nx-1
for iy = 2:ny-1
k = (iy-1)*nx + ix;
index = index + 1;
irc(index) = k;
icc(index) = k;
mat(index) = -dy/dx - dy/dx - dx/dy - dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k+1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k-1;
mat(index) = dy/dx;
index = index + 1;
irc(index) = k;
icc(index) = k+nx;
mat(index) = dx/dy;
index = index + 1;
irc(index) = k;
icc(index) = k-nx;
mat(index) = dx/dy;
%A(k,k) = -dy/dx - dy/dx - dx/dy - dx/dy;
%A(k,k+1) = dy/dx;
%A(k,k-1) = dy/dx;
%A(k,k+nx) = dx/dy;
%A(k,k-nx) = dx/dy;
b(k) = 0;
end
end
A = sparse(irc,icc,mat,nx*ny,nx*ny);
u = A\b;
%u
U = zeros(nx,ny);
for iy = 1:ny
for ix = 1:nx
k = (iy-1)*nx + ix;
U(ix,iy) = u(k);
end
end
[X,Y] = meshgrid(x,y);
surf(X, Y, U, 'EdgeColor','none');
end
function bc_value = bcfun(x,y)
if x==0 || y == 0
bc_value = 0;
return
end
if x==1
bc_value = y;
return
end
if y==1
bc_value = x;
return
end
end
##### 댓글 수: 2없음 표시없음 숨기기
skerdi hymeraj 2023년 9월 28일
sorry I have no idea how to use this site

댓글을 달려면 로그인하십시오.

### 카테고리

Help CenterFile Exchange에서 Matrix Indexing에 대해 자세히 알아보기

R2022a

### Community Treasure Hunt

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

Start Hunting!

Translated by