python 中 PDE 的复值 ODE

Question

我有一个来自 PDE 问题的复值系统，Python 中的 odeint() 无法处理它。我写了一个 RK4 模块来解决我的系统。它似乎有效，但是，计算值显然不正确。在第二个时间步，整个计算值为零。这是我的代码：

import matplotlib.pyplot as plt
import numpy as np
import drawnow
import time
import math

### Parameters ###
L       = 20
n       = 64
delta_t = 1.
tmax    = 10
miu     = 1e-6
x2      = np.linspace(-L/2,L/2, n+1)
x       = x2[:n]                         # periodic B.C. #0 = #n

kx1     = np.linspace(0,n/2-1,n/2)
kx2     = np.linspace(1,n/2,  n/2)
kx2     = -1*kx2[::-1]
kx      = (2.*math.pi/L)*np.concatenate((kx1,kx2)); kx[0] = 1e-6
ky      = kx;                                       y     = x

X, Y    = np.meshgrid(x, y)
KX,KY   = np.meshgrid(kx,ky)
K       = KX**2 + KY**2
K2      = np.reshape(K, n**2,1)

### Initial Condition ###
vorticity = np.exp(-0.25*X**2 - 2.*Y**2)
wt        = np.fft.fft2(vorticity)
wt2       = np.reshape(wt, n**2, 1)     # wt2 is initial condition
### Define ODE ###  
def yprime(t,rhs):  
    global miu, K, K2,n,KX, KY, wt2, wt 
    psit = -wt/ K   
    psix = np.real(np.fft.ifft2(1j*KX*psit))
    psiy = np.real(np.fft.ifft2(1j*KY*psit))
    wx   = np.real(np.fft.ifft2(1j*KX*wt))
    wy   = np.real(np.fft.ifft2(1j*KY*wt))
    rhs  = -miu * K2 * wt2 + np.reshape(np.fft.fft2(wx*psiy - wy*psix), n**2,1)
    return rhs

def RK4(domain,wt2,tmax):
    w = np.empty((tmax+1,n**2))
    w = w + 0j
    t = np.empty(tmax+1)        # length
    w[0,:] = wt2                # enter initial conditions in y
    t[0]   = domain[0]
    for i in range(1,tmax):
        t[i+1] = t[i]+delta_t
        w[i+1,:] = RK4Step(t[i], w[i,:],delta_t)
    return w

def RK4Step(t,w,delta_t):
    k1 = yprime(t,w)
    k2 = yprime(t+0.5*delta_t, w+0.5*k1*delta_t)
    k3 = yprime(t+0.5*delta_t, w+0.5*k2*delta_t)
    k4 = yprime(t+delta_t,     w+k3*delta_t)
    return w + (k1+2*k2+2*k3+k4)*delta_t/6. 

### Prediction ###
TimeStart = 0.
TimeEnd   = tmax+1
TimeSpan  = np.arange(TimeStart, TimeEnd, delta_t)
wt2_sol   = RK4(TimeSpan, wt2, tmax)

for i in TimeSpan:
    w = np.real(np.fft.ifft2(np.reshape(wt2_sol[i,:], (n, n))))
    plt.pcolor(X,Y,w,shading = 'interp',cmap='jet')
    drawnow
    time.sleep(0.2)
    plt.show()

知道为什么它不起作用吗？另外，我喜欢根据解决方案制作一个短视频。函数 'drawnow' 和 'time.sleep() 在这里似乎不起作用。

谢谢！

Answer 1

我的清理版本。更改内部步数不会改变结果的质量。

使 Runge-Kutta 求解器（更）通用，输入时间数组 (times[0],y0) 作为 IVP 的初始点
将def yprime(t,rhs):替换为def yprime(t,wt):，因为wt是你的状态变量，rhs是结果。所以wt是yprime中的局部变量。在return语句中直接汇编消除rhs。
删除所有reshape操作，作用于二维数组的向量space，numpy擅长将矩阵视为其他类型的向量
为预生成的图像序列添加matplotlib.animate。教程似乎比基于函数的动画更简单
在 kx 的生成中使用 arange 替换 linspace。更好的选择可能是使用 fftshift 交换频率数组的一半

.

import numpy as np
import math
from matplotlib import pyplot as plt
from matplotlib import animation

#----- Numerical integration of ODE via fixed-step classical Runge-Kutta -----

def RK4Step(yprime,t,y,dt):
    k1 = yprime(t       , y          )
    k2 = yprime(t+0.5*dt, y+0.5*k1*dt)
    k3 = yprime(t+0.5*dt, y+0.5*k2*dt)
    k4 = yprime(t+    dt, y+    k3*dt)
    return y + (k1+2*k2+2*k3+k4)*dt/6. 

def RK4(yprime,times,y0):
    y = np.empty(times.shape+y0.shape,dtype=y0.dtype)
    y[0,:] = y0                 # enter initial conditions in y
    steps = 4
    for i in range(times.size-1):
        dt = (times[i+1]-times[i])/steps
        y[i+1,:] = y[i,:]
        for k in range(steps):
            y[i+1,:] = RK4Step(yprime, times[i]+k*dt, y[i+1,:], dt)
    return y

#====================================================================
#----- Parameters for PDE -----
L       = 20
n       = 64
delta_t = 1.
tmax    = 10
miu     = 1e-6

#----- Constructing the grid and kernel functions
x2      = np.linspace(-L/2,L/2, n+1)
x       = x2[:n]                         # periodic B.C. #0 = #n
y       = x

kx     = np.linspace( -n/2  , n/2-1, n)
kx      = (2.*math.pi/L)*np.concatenate((np.arange(0,n/2),np.arange(-n/2,0))); 
kx[0] = 1e-6
ky      = kx;                                       

X, Y    = np.meshgrid(x, y)
KX,KY   = np.meshgrid(kx,ky)
K       = KX**2 + KY**2

#----- Initial Condition -----
vorticity = np.exp(-0.25*X**2 - 2.*Y**2)
wt0       = np.fft.fft2(vorticity)

#----- Define ODE -----
def wprime(t,wt):  
    global miu, K, K2,n,KX, KY 
    psit = -wt / K   
    psix = np.real(np.fft.ifft2(1j*KX*psit))
    psiy = np.real(np.fft.ifft2(1j*KY*psit))
    wx   = np.real(np.fft.ifft2(1j*KX*wt))
    wy   = np.real(np.fft.ifft2(1j*KY*wt))
    return  -miu * K * wt + np.fft.fft2(wx*psiy - wy*psix)

#====================================================================
#----- Compute the numerical solution -----
TimeStart = 0.
TimeEnd   = tmax+delta_t
TimeSpan  = np.arange(TimeStart, TimeEnd, delta_t)
wt_sol    = RK4(wprime,TimeSpan, wt0)

#----- Animate the numerical solution -----
fig = plt.figure()
ims = []
for i in TimeSpan:
    w = np.real(np.fft.ifft2(wt_sol[i,:]))
    im = plt.pcolor(X,Y,w,edgecolors='none',cmap='jet')
    ims.append([im])

ani = animation.ArtistAnimation(fig, ims, interval=50, blit=True,
    repeat_delay=1000)

#ani.save('PDE-animation.mp4')

plt.show()

python 中 PDE 的复值 ODE

Complex-valued ODE from PDE in python

python

ode