I'm trying to solve a coupled non-linear system of 4 PDEs, here's the code:
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
However it gives the following error and does not really move from the initial point:
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
And the final solution G is an array of 3996 size, please help me understand this. Thanks!
