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

Answer 1

[Torsten]

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); 

[Tristofani Agasta]

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]

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]

Thank you for your explanation! I think I can understand how the iterations work now

[Torsten]

Which iterations do you mean ? The loops to set up A and B ?

[Tristofani Agasta]

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

Answer 2

[Torsten]

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

[skerdi hymeraj]

sorry I have no idea how to use this site

[Torsten]

The internet doesn't forget: