2016-03-07 1 views
2

Это модифицированная задача из книги численных вычислений-Кинкейда, глава 15. (не физика)конечных разностей для 2D Parabolic PDE

Как правильно реализовать граничные условия? Условия:

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

Я не делаю это правильно, похоже. Мой код ниже.

Я пытаюсь написать код Fortran для решения уравнения теплоты (параболического) 2D, используя Finite Differences. Когда я распечатываю свои результаты, я получаю разные результаты и «NaN». Кажется, я не правильно определяю граничные условия. Я правильно сделал код в 1 измерении, пытаясь обобщить его на две части. У меня проблемы на границе.

Примечание: i,j предназначены для положений x и y соответственно, т. Е. Для цикла времени. nx,ny,W - количество точек сетки в направлении x, y и времени соответственно. Lx,Ly и tmax - это размер положения и временных интервалов для сетки. Этап положения (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 
+0

Здесь мы действительно не решаем проблемы физики, извините. Если у вас есть конкретная проблема с кодированием, мы, вероятно, сможем помочь, но этот тип проблемы, вероятно, выходит за рамки для stackoverflow. – Ross

+0

Сказав это, я могу предложить пару вопросов, которые могут помочь вам на вашем пути. Вы уверены, что правильно устанавливаете граничные условия? Похоже, вы только устанавливаете u (0,0) и u (nx, ny), но затем, например, используете u (0,1: ny-1). – Ross

+0

Кроме того, поскольку похоже, что вы пытаетесь учиться, я рекомендую вам отступы вашего кода. Очень трудно прочитать неверно отформатированный код, а отступом для циклов и условностей является ключевая часть этого. – Ross

ответ

2

Как Росс отмечает, вы не в полной мере определены граничные условия для ребер 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.

Внутренние пункты перезаписываются позже.