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

我的清理版本。改變內部步驟的數量不會改變結果的質量。

• 充分利用龍格 - 庫塔求解器（更多）通用，輸入時間陣列(times[0],y0)是的IVP

• 初始點與def yprime(t,wt):替換def yprime(t,rhs):，由於wt是你的狀態變量，並rhs是結果。所以wt是yprime中的局部變量。在return語句中通過直接組裝消除rhs。

• 刪除所有reshape操作中，二維陣列的矢量空間的行爲，numpy的是在治療矩陣作爲一些其他類型的矢量

• 的添加matplotlib.animate用於將所生成的預圖像序列好。這個教程似乎比基於功能的動畫更容易

• 與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()