二维抛物线 PDE 的有限差分
Finite differences for 2D Parabolic PDE
这是数值计算-Kincaid 的书第 15 章的一个修改后的问题。(不是物理学)
如何正确实施边界条件?条件是
u(0,y,t) = u(x,0,t) = u(nx,y,t) = u(x,ny,t) = 0.
看来我做的不对。我的代码如下。
我正在尝试编写 Fortran 代码来使用有限差分求解二维热(抛物线)方程。当我打印出我的结果时,我得到了不同的结果 'NaN'。看来我没有正确定义边界条件。我在一维中正确地编写了代码,试图将其归纳为两个我在边界上遇到了麻烦。
注意,i,j
分别是x和y位置do循环,m是时间do循环。 nx,ny,W
分别为x、y方向和时间上的网格点数。 Lx,Ly
和 tmax
是网格的位置和时间间隔的大小。 position(x,y) 步长和时间步长分别由 hx,hy,k
给出,下面的示例中 hx
和 hy
是相等的。我将我的解决方案存储在变量 u 和 v 中,如下所示。
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
正如 Ross 指出的那样,您还没有完全指定边 i=j=0
和 i=nx
以及 j=nx
的边界条件。只指定了您域的角落。
改变
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)
到
u(0,:)=0.0
u(nx,:)=0.0
u(:,0)=0.0
u(:,ny)=0.0
甚至
u=0.0
.
稍后会覆盖内部点。
这是数值计算-Kincaid 的书第 15 章的一个修改后的问题。(不是物理学)
如何正确实施边界条件?条件是
u(0,y,t) = u(x,0,t) = u(nx,y,t) = u(x,ny,t) = 0.
看来我做的不对。我的代码如下。
我正在尝试编写 Fortran 代码来使用有限差分求解二维热(抛物线)方程。当我打印出我的结果时,我得到了不同的结果 'NaN'。看来我没有正确定义边界条件。我在一维中正确地编写了代码,试图将其归纳为两个我在边界上遇到了麻烦。
注意,i,j
分别是x和y位置do循环,m是时间do循环。 nx,ny,W
分别为x、y方向和时间上的网格点数。 Lx,Ly
和 tmax
是网格的位置和时间间隔的大小。 position(x,y) 步长和时间步长分别由 hx,hy,k
给出,下面的示例中 hx
和 hy
是相等的。我将我的解决方案存储在变量 u 和 v 中,如下所示。
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
正如 Ross 指出的那样,您还没有完全指定边 i=j=0
和 i=nx
以及 j=nx
的边界条件。只指定了您域的角落。
改变
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)
到
u(0,:)=0.0
u(nx,:)=0.0
u(:,0)=0.0
u(:,ny)=0.0
甚至
u=0.0
.
稍后会覆盖内部点。