Trying to solve 2 dimensional Partial differential equation using Finite Difference Method
显示 更早的评论
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 个）
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 个评论
Tristofani Agasta
2022-3-22
Thank you for your correction. If I may ask, I tried to make 2d meshgrid with same code, is this the correct one?
% Showing Result
[X, Y] = meshgrid(x, y);
surface(X, Y, U); colormap
shading interp; axis ('equal')
xlim([0 max(x)]);ylim([0 max(y)])
xlabel('Sumbu X'); ylabel('Sumbu Y')
Torsten
2022-3-22
If you include a color map and choose more colors, you might see better.
The surface plot I used looked correct:
surf(X,Y,U)
I must admit that my experience in graphical representation is limited.
Tristofani Agasta
2022-3-22
Tristofani Agasta
2022-3-23
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
类别
在 帮助中心 和 File Exchange 中查找有关 Vector Volume Data 的更多信息
2022-3-21
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
