周期边界条件-有限差分
周期边界条件-有限差分
提问于 2016-03-09 22:29:35
嗨，下面我有一个解决非线性耦合PDE问题的代码，但是我需要实现周期性边界条件。周期性边界条件困扰着我，我应该在我的代码中添加什么来执行周期性边界条件？根据下面的模块算法建议进行更新。
注意，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谢谢你的阅读，如果我应该以更恰当的方式重新编排我的问题，请告诉我。
回答 1
Stack Overflow用户
回答已采纳
发布于 2016-03-10 10:43:54
PDE离散化中边界条件的一个选择是使用鬼(晕)单元(网格点)。对于周期BC，它可能不是最聪明的，但它可以用于所有其他边界条件类型。
因此，您将数组声明为
real, dimension(-1:n) :: u,v,w,f,g,d但是你只在点0..n-1 (点n和0点相同)中求解PDE。您还可以从1.n中进行操作，并从0..n+1中声明数组。
然后你设置
u(-1) = u(n-1)和
u(n) = u(0)所有其他数组也是如此。
在每一时间步骤中，您都会再次为u和f或在解决方案期间修改的所有其他字段设置这个值：
do m=1,A
u(-1) = u(n-1)
u(n) = u(0)
f(-1) = f(n-1)
f(n) = f(0)
do i=0,n-1 !Discretization equation for all times after the 1st step
u(i)=...
f(i)=...
end do
end do所有这些假设都是显式时间离散化和有限差分空间离散化，并且假设x(0) = 0和x(n) = 1是你的边界点。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/35903982
复制相关文章
相似问题
腾讯云开发者
