I am totally new to coding and I want to solve these 5 differential equations numerically. I took a python template and applied it to my case. Here's the simplified version of what I wrote:
import numpy as np
from math import *
from matplotlib import rc, font_manager
import matplotlib.pyplot as plt
from scipy.integrate import odeint
#Constants and parameters
alpha=1/137.
k=1.e-9
T=40.
V= 6.e-6
r = 6.9673e12
u = 1.51856e7
#defining dy/dt's
def f(y, t):
A = y[0]
B = y[1]
C = y[2]
D = y[3]
E = y[4]
# the model equations
f0 = 1.519e21*(-2*k/T*(k - (alpha/pi)*(B+V))*A)
f1 = (3*B**2 + 3*C**2 + 6*B*C + 2*pi**2*B*T + pi**2*T**2)**-1*(-f0*alpha/(3*pi**3) - 2*r*(B**3 + 3*B*C**2 + pi**2*T**2*B) - u*(D**3 - E**3))
f2 = u*(D**3 - E**3)/(3*C**2)
f3 = -u*(D**3 - E**3)/(3*D**2)
f4 = u*(D**3 - E**3)/(3*E**2) + r*(B**3 + 3*B*C**2 + pi**2*T**2*B)/(3*E**2)
return [f0, f1, f2, f3, f4]
# initial conditions
A0 = 2.e13
B0 = 0.
C0 = 50.
D0 = 50.
E0 = C0/2.
y0 = [A0, B0, C0, D0, E0] # initial condition vector
t = np.linspace(1e-15, 1e-10, 1000000) # time grid
# solve the DEs
soln = odeint(f, y0, t, mxstep = 5000)
A = soln[:, 0]
B = soln[:, 1]
C = soln[:, 2]
D = soln[:, 3]
E = soln[:, 4]
y2 = [A[-1], B[-1], C[-1], D[-1], E[-1]]
t2 = np.linspace(1.e-10, 1.e-5, 1000000)
soln2 = odeint(f, y2, t2, mxstep = 5000)
A2 = soln2[:, 0]
B2 = soln2[:, 1]
C2 = soln2[:, 2]
D2 = soln2[:, 3]
E2 = soln2[:, 4]
y3 = [A2[-1], B2[-1], C2[-1], D2[-1], E2[-1]]
t3 = np.linspace(1.e-5, 1e1, 1000000)
soln3 = odeint(f, y3, t3)
A3 = soln3[:, 0]
B3 = soln3[:, 1]
C3 = soln3[:, 2]
D3 = soln3[:, 3]
E3 = soln3[:, 4]
#Plot
rc('text', usetex=True)
plt.subplot(2, 1, 1)
plt.semilogx(t, B, 'k')
plt.semilogx(t2, B2, 'k')
plt.semilogx(t3, B3, 'k')
plt.subplot(2, 1, 2)
plt.loglog(t, A, 'k')
plt.loglog(t2, A2, 'k')
plt.loglog(t3, A3, 'k')
plt.show()
I get the following error:
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 = 0.2135341098625E-06 R2 = 0.1236845248713E-22
For some set of parameters, upon playing with mxstep
in odeint (also tried hmin
and hmax
but didn't notice any difference), although the error persists, my graphs look good and are not impacted, but most of the times they are.
Sometimes the error I get asks me to run with the odeint option full_output=1
and in doing so I obtain:
A2 = soln2[:, 0]
TypeError: tuple indices must be integers, not tuple
I don't get what this means when searching for it.
I would like to understand where the problem lies and how to solve it. Is odeint even suitable for what I'm trying to do?