【发布时间】：2016-06-24 13:25:16
【问题描述】：
您好，我有一个解决非线性耦合偏微分方程的代码。但是我需要实现周期性边界条件。周期性边界条件困扰着我，我应该在代码中添加什么来强制执行周期性边界条件？根据以下模块化算术建议进行了更新。
注意，t>=0 且 x 在区间 [0,1] 内。这是耦合方程，下面我提供我的代码
其中 a, b > 0。
这些是初始条件，但现在我需要施加周期性边界条件。这些在数学上可以写成 u(0,t)=u(1,t) 和 du(0,t)/dx=du(1,t)/dx，同样适用于 f(x,t)。我对边界条件的 du/dx 实际上是偏导数。
我的代码在下面
program coupledPDE
integer, parameter :: n = 10, A = 20
real, parameter :: h = 0.1, k = 0.05
real, dimension(0:n-1) :: u,v,w,f,g,d
integer:: i,m
real:: t, R, x,c1,c2,c3,aa,b
R=(k/h)**2.
aa=2.0
b=1.0
c1=(2.+aa*k**2.-2.0*R)/(1+k/2.)
c2=R/(1.+k/2.)
c3=(1.0-k/2.)/(1.0+k/2.)
c4=b*k**2./(1+k/2.)
do i = 0,n-1 !loop over all space points except 0 and n
x = real(i)*h
w(i) = z(x) !u(x,0)
d(i) = z(x) !f(x,0)
end do
do i=0,n-1
ip=mod(i+1,n)
il=modulo(i-1,n)
v(i) = (c1/(1.+c3))*w(i) + (c2/(1.+c3))*(w(ip)+w(il)) -(c4/(1.+c3))*w(i)*((w(i))**2.+(d(i))**2.) !\partial_t u(x,0)=0
g(i) = (c1/(1.+c3))*d(i) + (c2/(1.+c3))*(d(ip)+d(il)) -(c4/(1.+c3))*d(i)*((w(i))**2.+(d(i))**2.) !\partial_t f(x,0)=0
end do
do m=1,A
do i=0,n-1
ip=mod(i+1,n)
il=modulo(i-1,n)
u(i)=c1*v(i)+c2*(v(ip)+v(il))-c3*w(i)-c4*v(i)*((v(i))**2.+(g(i))**2.)
f(i)=c1*g(i)+c2*(g(ip)+g(il))-c3*d(i)-c4*g(i)*((v(i))**2.+(g(i))**2.)
end do
print*, "the values of u(x,t+k) for all m=",m
print "(//3x,i5,//(3(3x,e22.14)))",m,u
do i=0,n-1
w(i)=v(i)
v(i)=u(i)
d(i)=g(i)
t=real(m)*k
x=real(i)*k
end do
end do
end program coupledPDE
function z(x)
real, intent(in) :: x
real :: pi
pi=4.0*atan(1.0)
z = sin(pi*x)
end function z
感谢阅读，如果我应该以更正确的方式重新格式化我的问题，请告诉我。
【问题讨论】：
-
只需在每个时间步中复制从 n 到 0 和从 1 到 n+1 的值。就是这样。
-
当然，只在网格
0,h,2h,...,1-h上声明所有内容，然后使用一些简单的模运算来移动索引会更干净。周期性边界被自动处理。更好的是使用整个数组和CSHIFT函数。 -
@VladimirF 感谢您的帮助，但我不明白您的意思。我在哪里复制这些值，你的意思是什么值？当我开始我的时间循环时，我会添加代码吗？ （我知道这是一个非常简单的实现想法，对于我的困惑感到抱歉）。再次感谢！
-
@RussF 也感谢您的回复。我如何在网格上声明它？我在掌握如何做到这一点时遇到了真正的问题。感谢您的帮助
-
只需将所有数组声明为
v[0:n-1]等。在循环中计算ip=mod(i+1,n)并将i+1的引用替换为ip。i-1类似。如果您将循环定义为do i = 0,n-1，一切都会自然处理。如果您确实需要 x=1 处的值，那么只需在i循环之后从 x=0 复制值即可。
标签： math fortran fortran90 numerical-methods pde