CODEWITHSUNDEEP
mathfortranfortran90gfortrandifferential-equations
Joe
Joe

Reputation: 33

Finite differences for 2D Parabolic PDE

This is a modified problem from Numerical Computing-Kincaid's book, chapter 15. (not physics)

How can I properly implement the boundary conditions? The conditions are 

u(0,y,t) = u(x,0,t) = u(nx,y,t) = u(x,ny,t) = 0.

I am not doing it correctly, it seems. My code is below.

I am trying to write a Fortran code to solve the 2D heat (parabolic) equation using Finite Differences. When I print out my results, I get divergent results and 'NaN'. It seems I am not defining the boundary conditions correctly. I properly did the code in 1 dimension, trying to generalize it in two I have troubles at the boundary.

Note, i,j are for the x and y position do loops respectively, m is for time do loop. nx,ny,W are the number of grid points in x , y direction and time respectively. Lx,Ly and tmax are the size of the position and time intervals for the mesh. The position(x,y) steps and time steps are given by hx,hy,k respectively, and hx and hy are equal for the example below. I store my solutions in the variables u and v as shown below.

program parabolic2D

implicit none

integer :: i,j,m 
integer, parameter :: nx=10., ny=10., W=21. 
real, parameter :: Lx=1.0, Ly=1.0, tmax=0.1 
real :: hx,hy,k,pi,pi2,R,t 
real, dimension (0:nx,0:ny) :: u,v 


hx=(Lx-0.0)/nx 
hy=(Ly-0.0)/ny 
k=(tmax-0.0)/W
R=k/hx**2.
u(0,0)=0.0; v(0,0)=0.0; u(nx,ny)=0.0; v(nx,ny)=0.0 !boundary conditions u(0,0,t)=0=u(nx,ny,t)
pi=4.0*atan(1.0) 
pi2=pi*pi



do i=1,nx-1
do j=1,ny-1
u(i,j)=sin(pi*real(i)*hx)*sin(pi*real(j)*hy)  !initial condition
end do
end do

do m=1,W

do i=1,nx-1
do j=1,ny-1
v(i,j) = R*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1))+(1-4*R)*u(i,j) !Discretization for u(x,y,t+k)
end do
end do

t = real(m)*k ! t refers to time in the problem.

do i=1,nx-1
do j=1,ny-1
u(i,j)=v(i,j) !redefining variables.
end do
end do
write(*,*) 'for all times m, this prints out u(x,y,t)',m,((u(i,j),i=0,nx),j=0,ny)

end do

end program parabolic2D

Upvotes: 3

Views: 545

Answers (1)

RussF
RussF

Reputation: 721

As Ross points out, you haven't fully specified the boundary conditions for the edges i=j=0 and i=nx and j=nx. Only the corners of your domain have bee specified.

Change

u(0,0)=0.0; v(0,0)=0.0; u(nx,ny)=0.0; v(nx,ny)=0.0 !boundary conditions u(0,0,t)=0=u(nx,ny,t)

to

u(0,:)=0.0 u(nx,:)=0.0 u(:,0)=0.0 u(:,ny)=0.0

or even

u=0.0.

The interior points are overwritten later.

Upvotes: 3

Related Questions