CODEWITHSUNDEEP
pythonpython-3.xodenumerical-integrationodeint
TQM
TQM

Reputation: 89

Scipy's Odeint Returns Zeros After Initial Time Step Solutions

In attempting to numerically solve for a system of ODEs, for some reason in the 3rd step of the solution, it turns to all zeros. I would think this implies an unstable solution, but in general I thought that case would return infinities or nan, not zeros. Why is this occurring, and how can I fix it? I am getting the following warning on running odeint (though its still completing): 

Warning

python\python37\lib\site-packages\scipy\integrate\odepack.py:236: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. 


 warnings.warn(warning_msg, ODEintWarning)

Code

import scipy as sp
from scipy.integrate import odeint
import numpy as np

def ODEEq(rho, t, D, B, Nd):
    return (D*np.gradient(np.gradient(rho))-B*(Nd*rho+(rho**2)))
def Solver(delz, delt, L, T, G0, D, B, Nd, Sf, Sb, unitMaker=1e-9):
    xm1Ghost = G0[1] - (Sf*delz*G0[0]/D)
    xLp1Ghost = G0[-2] - (Sb*delz*G0[-1]/D)
    rho = G0
    np.insert(rho, 0, xm1Ghost)
    np.append(rho, xLp1Ghost)

    sol = odeint(ODEEq, rho, tspace, args=(D, B, Nd), full_output=1)
    return sol

L = 4000
timeRange = 2000
D = 1.639e18
B = 2e13
Nd = 0
delz = L/(10*L)
delt = timeRange/(10*timeRange)
G0 = 640*np.exp((-6.4e-2)*np.linspace(0,1000,L))
Sf = 0.64
Sb = 0.64

sol = Solver(delz, delt, L, timeRange, G0, D, B, Nd, Sf, Sb)

Result (sol):

(array([[6.40000000e+02, 6.29838965e+02, 6.19839253e+02, ...,
    1.05982469e-25, 1.04299825e-25, 1.02643897e-25],
   [3.64489112e+01, 3.64037414e+01, 3.63580081e+01, ...,
    1.64712980e-25, 1.61424855e-25, 1.58213241e-25],
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
   ...,
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]), {'hu': array([2.24035532e-005, 8.91443222e-312, 4.78098028e+004, ...,
   3.21768899e-021, 3.11632803e-021, 3.01816006e-021]), 'tcur': array([1.28018934e-003, 8.91443222e-312, 2.65545333e-021, ...,
   4.31714361e-049, 4.15062501e-049, 3.98753839e-049]), 'tolsf': array([ 0.00000000e+000,  8.91443222e-312, -4.85590765e+004, ...,
    5.32418663e-009,  5.23965656e-009,  5.15646853e-009]), 'tsw': array([2.54756216e-004, 8.91443222e-312, 4.83671586e-009, ...,
   5.52645557e-023, 5.40880326e-023, 5.29122707e-023]), 'nst': array([        500,         420,  2051404512, ...,  -901543794,
    -382763050, -1396204354], dtype=int32), 'nfe': array([     52989,        420, 2051404512, ...,  896695326, -858716839,
   -564525364], dtype=int32), 'nje': array([        13,        420,  604126608, ..., 1834441325, 1211124015,
   1163472949], dtype=int32), 'nqu': array([         5,        420,  414912656, ..., -441395200, 1046325818,
    744409088], dtype=int32), 'imxer': -1, 'lenrw': 16036022, 'leniw': 4020, 'mused': array([         2,        420,  414912656, ...,   67999744, 1022034837,
    973886464], dtype=int32), 'message': 'Excess work done on this call (perhaps wrong Dfun type).'})

Upvotes: 0

Views: 337

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25972

Your system appears to be u_t = D*u_xx - B*(Nd*u+u**2) with boundary conditions D*u_x(0)=Sf*u(0)/2 and D*u_x(L)=Sb*u(L)/2, which is a diffusion process with a non-linear local dynamic that has a stable equilibrium at u=0 and an unstable equilibrium at u=-Nd.

  • As you have free or at least unspecified boundaries per your implementation, it is no enigma that the system converges to its stable equilibrium everywhere.

  • Your setup is incomprehensible, what is the length of the space domain, what its subdivision? You set delz = 0.1, but the spacing in np.linspace(0,1000,L) is 1000/3999, about 0.25

  • You do not apply the space segmentation to the gradient computation, you would either need to pass delz as argument or divide by its square afterwards.

  • You need to apply the boundary condition, resp. the computation of the ghost values, to every evaluation of the ODE function.

To deeper discuss how to translate your PDE into a method-of-lines system, please post equations and your ideas so far in scientific computing SE, as that mathematical content is off-topic here. As far as I can see, there is no programming error involved, only the mathematical-numerical characteristics of this, likely wrongly constructed, ODE system.

Upvotes: 1

Related Questions