Computational Science Asked by b-dog on June 13, 2021
I am trying out odeint
and received the error
‘Excess work done on this call (perhaps wrong Dfun type).’.
The values returned are also super sensitive to small changes of the initial value. Some initial values don’t run into this problem whereas others do.
rho = 0.1
beta = 0.01
N = 50
gamma = 0.20
theta = 0.01
delta = 0.1
alpha = 1000
S = 1000
W = 180
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def model(L,t):
R = (L * delta * alpha * 20) / (2* theta*S*N*W)
dLdt= rho*L-beta*N+(L*delta*alpha*R*W*20)/(N*S)
return dLdt
L0 = 20
# time points
t = np.linspace(0, 500, 1000)
# solve ODE
L, infodict = odeint(model,L0,t,full_output=True)
infodict['message']
The output L
goes from 20
to 100
and then spikes to 1.62e+011
before going to zero again.
Essentially, by combining all the constant factors, your ODE is $$ frac{dL}{dt} = -A + B L + C L^2 $$ With the initial $L(0)=L_0$ large enough, the positive feed-back of the quadratic term drives the equation towards a dynamic blow-up in finite time, as you observed with the values getting very large. This is a property of the exact solution, so it is no surprise that an accurate numerical solution will replicate this.
Your equation will only give bounded solutions if the initial point is below the positive root of the quadratic polynomial on the right side.
You have
A = beta*N
B = rho
C = ((delta * alpha * 20)/(N*S))**2/(2*theta)
print([A,B,C])
print(np.roots([C,B,-A]))
giving coefficients [0.5, 0.1, 0.08]
and roots of the quadratic [-3.20194102 1.95194102]
. For values slightly larger than that you will still reach the end of the integration interval before the blow-up point, however the values will be getting larger.
Correct answer by Lutz Lehmann on June 13, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP