滚转“变函数”为循环的等效程序及经典方法

Question

我有一个基本循环(实际上是时间循环，其中更新数组的值)：

for i in range(1,nt):
    #Using roll
    u = u - cfl/2*(roll(u,-1)- roll(u,1))
    # Update time
    t = t+dt

利用该回路和roll函数，仿真运行良好。

我精确地说我必须有周期性的边界条件(x[0] = x[n-1])

现在，我试图得到相同的等价代码，但使用经典方法，我的意思是编写roll移位函数的显式代码。

根据先前在另一个论坛上的一篇文章，有人建议这样做：

for i in range(1,nt):
   utemp1 = u[0] - cfl/2*(u[nx-1] - u[1])
   utemp2 = u[nx-1] - cfl/2*(u[nx-2] - u[0])
   u[1:nx-2] = u[1:nx-2] - cfl/2*(u[0:nx-3] - u[2:nx-1])
   u[0] = utemp1
   u[nx-1] = utemp2

但是，这个“等价代码”不会产生与第一个版本(使用roll python函数)相同的结果。

我想了解我上面的作业有什么不对。这似乎是一个覆盖的值或边界上的错误更新的问题。

更新1

以下代码对我来说似乎是正确的：

for i in range(1,nt):
   utemp1=u[0] - cfl/2*(u[nx-1] - u[1])
   utemp2=u[nx-1] - cfl/2*(u[nx-2] - u[0])
   u[1:nx-2] = u[1:nx-2] - cfl/2*(u[0:nx-3] - u[2:nx-1])
   u[0] = utemp1
   u[nx-1] = utemp2

除了我不确定是否应该倒置最后两项任务(由于定期边界条件)以外：

   u[0] = utemp1
   u[nx-1] = utemp2

这个"roll“Python函数的”等效代码“必须”正常“工作，这只是而且不幸的是我的印象，因为它没有产生与简单解决方案相同的结果：

for i in range(1,nt):
    #Using roll
    u = u - cfl/2*(roll(u,-1)- roll(u,1))

为什么我的显式解决方案( "roll“python函数的显式工作)不能工作？

更新2

根据@Engineero的请求，我提供了一个可以简单使用Python运行的小代码：

#
# Equation d'advection: u,t + V*u,x = 0
# centre FTCS (Forward Time Centered Space)
#

import matplotlib.pyplot as plt
from pylab import *

# Speed
V = 1  
L = 1
# analytical solution --------------------------
def uexacte(t,x):
    return sin(2*pi*(x-V*t)/L)

# 1. Centre FTCS (Forward Time Centered Space)
cfl = 0.1
nx = 128
tend = 1

#
dx = L/(nx-1.)
dt = cfl*dx/V
nt = int(tend/dt)+1
print "CFL=%5.2f tend=%4.1f --> %i iterations en temps"%(cfl,tend,nt)

# Arrays
x = linspace(0,L,nx)
# Bounadry condition
u0 = uexacte(0,x)

# Starting solution
t=0.0 ; u=copy(u0)

# Time loop
for i in range(1,nt):
    # FTCS

    # Using classical approach
    utemp1=u[0] - cfl/2*(u[nx-1] - u[1])
    utemp2=u[nx-1] - cfl/2*(u[nx-2] - u[0])
    u[1:nx-1] = u[1:nx-1] - cfl/2*(u[0:nx-2] - u[2:nx])
    u[0] = utemp1
    u[nx-1] = utemp2
 
    #Using roll
    #u = u - cfl/2*(roll(u,-1)- roll(u,1))
    
    # Update time
    t = t+dt

    # Plot numerical solution and analytical
    #if (i % 2 == 0):
    plt.gcf().clear()
    # Numerical
    plt.ylim(ymin=-1.1, ymax=1.1)
    plt.title('FTCS Scheme - Unconditionally unstable - CFL=%5.2f'%cfl)
    plt.xlabel('x')
    plt.ylabel('u')
    plt.plot(x,u,'-o',color='r',label='Numerical')
    # Analytical
    uexact= uexacte(t,x)
    plt.plot(x,uexact,'--',label='Analytical')
    plt.draw()
    plt.legend()
    plt.pause(0.1)

plt.pause(30)

如果我使用我的经典方法(“roll”的显式工作)，我就得不到解析解。

下面是动画过程中捕获的一个示例：

​

​

在代码开始时，uexacte函数给出了解析解：

sin(2*pi*(x-V*t)/L)

另一方面，通过使用"roll“函数，在解析解和数值解之间的动画过程中得到了相同的结果(直到某个点)。

更新3

我一直试图通过以下方式比较这两个公式(与roll )的关系，使之成为roll的经典方法：

# Time loop
for i in range(1,nt):
    # FTCS

    # Using classical approach
    utemp1 = u[0] - cfl/2*(u[nx-1] - u[1])
    utemp2 = u[nx-1] - cfl/2*(u[nx-2] - u[0])
    u[1:nx-1] = u[1:nx-1] - cfl/2*(u[0:nx-2] - u[2:nx])
    u[0] = utemp1
    u[nx-1] = utemp2
    print "array u"
    print(u)

    #Using roll
    utest = utest - cfl/2*(roll(utest,-1) - roll(utest,1))
    print "array utest"
    print(utest)

在执行死刑的时候，我得到：

array u
[ 0.04330127  0.82272413 -0.90932667  0.04330127]
array utest
[-0.04330127  0.90932667 -0.82272413 -0.04330127]
array u
[ 0.08227241  0.77509274 -0.94829782  0.09093267]
array utest
[-0.09093267  0.94829782 -0.77509274 -0.08227241]
array u
[ 0.11648042  0.72356422 -0.98250582  0.14246118]
array utest
[-0.14246118  0.98250582 -0.72356422 -0.11648042]

我刚开始第一次迭代，两个数组的值就不相对应了，我不明白为什么.

更新4

好的，我想我已经找出了这个问题：

如果我在一次操作中执行“多索引”的赋值(我指一次)：

u[1:nx-2] = u[1:nx-2] - cfl/2.0*(u[2:nx-1] - u[0:nx-3])

赋值执行得很差，阵列u的值计算也不太好。

如果我用以下更简单的方法测试“一次和一次分配”(只有两个术语)：

u[1:nx-2] = u[1:nx-2] - u[2:nx-1]

正确计算了阵列1:nx-2的u值。

如何执行这种自然的，当然过于简单的操作形式：

u[1:nx-2] = u[1:nx-2] - cfl/2.0*(u[2:nx-1] - u[0:nx-3])

?？?(也许这个因素不适用于不同的u[2:nx-1] - u[0:nx-3]，并且在u1:nx-2中也有这种差异，这可能意味着问题)。

Answer 1

两个同等职能：

nt=7
nx=10
cfl=1.1
def old(u0,cfl=-cfl):
    u=u0.copy()
    for i in range(1,nt):
        utemp1=u[0] - cfl/2*(u[nx-1] - u[1])
        utemp2=u[nx-1] - cfl/2*(u[nx-2] - u[0])
        u[1:nx-1] = u[1:nx-1] - cfl/2*(u[0:nx-2] - u[2:nx])
        u[0] = utemp1
        u[nx-1] = utemp2
    return u


def new(u0):
    u=u0.copy().astype(float)        
    for i in range(1,nt):
        u -=  cfl/2*(roll(u,-1)- roll(u,+1))
    return u

指数、滚动标志和v类型的错误:它必须是浮动的。

运行：

In [521]: u0
Out[521]: array([ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.])

In [522]: old(u0)
Out[522]: 
array([-15.06651094,  22.56142656,  28.93808047,  11.76161172,
         0.41970625,   1.41970625,  -9.92219922,   9.25426953,
        15.63092344, -19.99701406])

In [523]: new(u0)
Out[523]: 
array([-15.06651094,  22.56142656,  28.93808047,  11.76161172,
         0.41970625,   1.41970625,  -9.92219922,   9.25426953,
        15.63092344, -19.99701406])