瞭解scipy.integrate，odeint包裝odeintw

我試圖解決的4個偏微分方程耦合非線性系統，下面的代碼：瞭解scipy.integrate，odeint包裝odeintw

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.integrate import odeint 
from odeintw import odeintw 

# constants 
global rH, L, q, p, Q, mu 
rH = 1.0 ; L = 1.0 ; p = 2.0 ; q = 1.0 ; Q = 1.71 ; mu = Q*rH/L**2 

def dfunc(G, kx, ky, w, r): 

try: 
    z = 1.0/(rH - r) 

except ZeroDivisionError: 
print "Trying to divide zero" 

phi = mu*(1.0 - rH/r) 
fr = 1.0 - (rH/r)**3 
scale = L**2 * np.sqrt(fr)/r**2 + 0j 

vp = (w + q*phi)/(np.sqrt(fr)) + p*mu*rH/r + 0j 
vn = (w + q*phi)/(np.sqrt(fr)) - p*mu*rH/r + 0j 

# G11 -> G1, G12 -> G2, G21 -> G3, G22 -> G4 

G1, G2, G3, G4 = G 

G11 = (vp - kx + (vn + kx)*G1**2 - ky*G1*G2 - ky*G1*G3 
+ (vp-kx)*G2*G3) 
G12 = (ky + (vn + kx)*G1*G2 - ky*G2**2 - ky*G1*G4 
+ (vp-kx)*G2*G4) 
G21 = (ky + (vn + kx)*G1*G3 - ky*G1*G4 - ky*G3**2 
+ (vp-kx)*G3*G4) 
G22 = (vn + kx + (vn + kx)*G2*G3 - ky*G2*G4 - ky*G2*G3 
+ (vp-kx)*G4**2) 

return np.multiply(scale, np.array([ G11, G12, G21, G22 ])) 

# jaccobian 

def jac(G, kx, ky, w, r): 

G1, G2, G3, G4 = G 

jac = np.array([ 
    [2(vn+kx)*G1 - ky*(G2+G3), -ky*G1 + (vp-kx)*G3, -ky*G1+(vp-kx)*G2, 0 ], 
    [(vn+kx)*G2-ky*G4, (vn+kx)*G1-2*ky*G2+(vp-kx)*G4, 0, -ky*G1+(vp-kx)*G2 ], 
    [(vn+kx)*G3-ky*G4, 0, (vn+kx)*G1-2*ky*G3+(vp-kx)*G4, -ky*G1+(vp-kx)*G3 ], 
    [0, (vn+kx)*G3-ky*(G4+G3), (vn+kx)*G2-ky*G2, -ky*G2+2*(vp-kx)*G4 ] 
    ]) 

return np.multiply(scale, jac) 

# Initial value 
Gir = np.array([ 1.0j, 0, 0, 1.0j ]) 
r = np.arange(rH + 0.1, 1000*rH, 1) 

kx = 0 
ky = 0 
w = 0.0001 + 0.0000001j 

G, infodict = odeintw(dfunc, Gir, r, args=(kx,ky,w), Dfun=jac, full_output=True) 

print 'final', np.size(G), G

但是它提供了以下錯誤，並沒有真正從移動初始點：

lsoda-- warning..internal t (=r1) and h (=r2) are 
    such that in the machine, t + h = t on the next step 
    (h = step size). solver will continue anyway 
    In above, R1 = .1100000000000E+01 R2 = .3569471009368E-22

而最終的解決方案G是一個3996大小的數組，請幫我理解這一點。謝謝！

來源

2015-11-17 user40689

在函數dfunc和jac中，第二個參數必須是自變量。在你的情況下，自變量是r，這樣的定義應該開始

def dfunc(G, r, kx, ky, w):

和

def jac(G, r, kx, ky, w):

來源

2015-11-17 22:48:59

我現在覺得自己很蠢，謝謝！ – user40689

我可以請你在odeintw上看看這個問題嗎？ https://stackoverflow.com/questions/47295511/using-odeint-odeintw-for-complex-equations-when-initial-conditions-are-real – Joel

瞭解scipy.integrate，odeint包裝odeintw

回答

相關問題