should help you out I finally got a reasonable solution for the system of ODEs by providing the Jacobian matrix to the solver. See below for my working solution. import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d
# Data
time, qheat = np.loadtxt('timeq.txt', unpack=True)
# Calculate Temperatures
interp_qheat = interp1d(time, qheat)
def tc_dt(t, tc, ts, q):
"""
dTc/dt = (TsTc)/(Rc*Cc) + Q/Cc
"""
rc = 1.94
cc = 62.7
return ((ts  tc) / (rc * cc)) + q / cc
def ts_dt(t, tc, ts):
"""
dTs/dt = (TfTs)/(Ru*Cs)  (TsTc)/(Rc*Cs)
"""
rc = 1.94
ru = 3.08
cs = 4.5
tf = 298.15
return ((tf  ts) / (ru * cs))  ((ts  tc) / (rc * cs))
def jacobian(t, y):
"""
Given the following system of ODEs
dTc/dt = (TsTc)/(Rc*Cc) + Q/Cc
dTs/dt = (TfTs)/(Ru*Cs)  (TsTc)/(Rc*Cs)
determine the Jacobian matrix of the righthand side as
Jacobian matrix = [df1/dTc, df2/dTc]
[df1/dTs, df2/dTs]
where
f1 = (TsTc)/(Rc*Cc) + Q/Cc
f2 = (TfTs)/(Ru*Cs)  (TsTc)/(Rc*Cs)
"""
cc = 62.7
cs = 4.5
rc = 1.94
ru = 3.08
jc = np.array([
[1 / (cc * rc), 1 / (cs * rc)],
[1 / (cc * rc), 1 / (cs * ru)  1 / (cs * rc)]
])
return jc
def func(t, y):
"""
Righthand side of the system of ODEs.
"""
q = interp_qheat(t)
tcdt = tc_dt(t, y[0], y[1], q)
tsdt = ts_dt(t, y[0], y[1])
return tcdt, tsdt
t0 = time[0]
tf = time[1]
sol = solve_ivp(func, (t0, tf), (298.15, 298.15), method='BDF', t_eval=time, jac=jacobian)
# Plot
fig, ax = plt.subplots(tight_layout=True)
ax.plot(sol.t, sol.y[0], label='tc')
ax.plot(sol.t, sol.y[1], label='ts')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Temperature [K]')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
plt.show()
Boards Message : 
You Must Login
Or Sign Up
to Add Your Comments . 
Share :
Tags :
Solution from SciPy solve_ivp contains oscillations for a system of firstorder ODEs ,
Solution ,
from ,
SciPy ,
solve_ivp ,
contains ,
oscillations ,
system ,
firstorder ,
ODEs ,

Solving a system of odes (with changing constant!) using scipy.integrate.odeint?
Date : March 29 2020, 07:55 AM
I wish did fix the issue. Yes, this is possible. In the case where a is constant, I guess you called scipy.integrate.odeint(fun, u0, t, args) where fun is defined as in your question, u0 = [x0, y0, z0] is the initial condition, t is a sequence of time points for which to solve for the ODE and args = (a, b, c) are the extra arguments to pass to fun. In the case where a depends on time, you simply have to reconsider a as a function, for example (given a constant a0): def a(t):
return a0 * t
def fun(u, t, a, b, c):
x = u[0]
y = u[1]
z = u[2]
dx_dt = a(t) * x + y * z # A change on this line: a > a(t)
dy_dt = b * (y  z)
dz_dt =  x * y + c * y  z
return [dx_dt, dy_dt, dz_dt]
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import scipy.integrate
tmax = 10.0
def a(t):
if t < tmax / 2.0:
return ((tmax / 2.0)  t) / (tmax / 2.0)
else:
return 1.0
def func(x, t, a):
return  (x  a(t))
x0 = 0.8
t = np.linspace(0.0, tmax, 1000)
args = (a,)
y = sp.integrate.odeint(func, x0, t, args)
fig = plt.figure()
ax = fig.add_subplot(111)
h1, = ax.plot(t, y)
h2, = ax.plot(t, [a(s) for s in t])
ax.legend([h1, h2], ["y", "a"])
ax.set_xlabel("t")
ax.grid()
plt.show()

How to solve a system of ODEs with scipy.integrate.odeint with a timedependent variable
Tag : python , By : Bjørn Lyngwa
Date : March 29 2020, 07:55 AM
To fix the issue you can do As far as I understood from the comments below your question, you try to incorporate measured data that can be noisy. Rather than plugging the data in directly, you can use these data to fit your timecourses. Here I show the outcome for your variable S: # zombie apocalypse modeling
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy.integrate import odeint
# solve the system dy/dt = f(y, t)
def f(y, t, paras):
Si = y[0]
Zi = y[1]
Ri = y[2]
try:
P = paras['P'].value
d = paras['d'].value
B = paras['B'].value
G = paras['G'].value
A = paras['A'].value
except:
P, d, B, G, A = paras
# the model equations (see Munz et al. 2009)
f0 = P  B * Si * Zi  d * Si
f1 = B * Si * Zi + G * Ri  A * Si * Zi
f2 = d * Si + A * Si * Zi  G * Ri
return [f0, f1, f2]
def g(t, x0, paras):
"""
Solution to the ODE x'(t) = f(t,x,p) with initial condition x(0) = x0
"""
x = odeint(f, x0, t, args=(paras,))
return x
def residual(paras, t, data):
x0 = paras['S0'].value, paras['Z0'].value, paras['R0'].value
model = g(t, x0, paras)
s_model = model[:, 0]
return (s_model  data).ravel()
# just for reproducibility reasons
np.random.seed(1)
# initial conditions
S0 = 500. # initial population
Z0 = 0 # initial zombie population
R0 = 0 # initial death population
y0 = [S0, Z0, R0] # initial condition vector
t = np.linspace(0, 5., 100) # time grid
P = 12 # birth rate
d = 0.0001 # natural death percent (per day)
B = 0.0095 # transmission percent (per day)
G = 0.0001 # resurect percent (per day)
A = 0.0001 # destroy percent (per day)
# solve the DEs
soln = odeint(f, y0, t, args=((P, d, B, G, A), ))
S = soln[:, 0]
Z = soln[:, 1]
R = soln[:, 2]
# plot results
plt.figure()
plt.plot(t, S, label='Living')
plt.plot(t, Z, label='Zombies')
plt.xlabel('Days from outbreak')
plt.ylabel('Population')
plt.title('Zombie Apocalypse  No Init. Dead Pop.; No New Births.')
plt.legend(loc=0)
plt.show()
# generate fake data
S_real = S[0::8]
S_measured = S_real + np.random.randn(len(S_real)) * 100
t_measured = t[0::8]
plt.figure()
plt.plot(t_measured, S_real, 'o', color='g', label='real data')
# add some noise to your data to mimic measurement erros
plt.plot(t_measured, S_measured, 'o', color='b', label='noisy data')
# set parameters including bounds; you can also fix parameters (use vary=False)
params = Parameters()
params.add('S0', value=S0, min=490., max=510.)
params.add('Z0', value=Z0, vary=False)
params.add('R0', value=R0, vary=False)
params.add('P', value=10, min=8., max=12.)
params.add('d', value=0.0005, min=0.00001, max=0.005)
params.add('B', value=0.01, min=0.00001, max=0.01)
params.add('G', value=G, vary=False)
params.add('A', value=0.0005, min=0.00001, max=0.001)
# fit model
result = minimize(residual, params, args=(t_measured, S_measured), method='leastsq') # leastsq nelder
# check results of the fit
data_fitted = g(t, y0, result.params)
plt.plot(t, data_fitted[:, 0], '', linewidth=2, color='red', label='fitted data')
plt.legend()
# display fitted statistics
report_fit(result)
plt.show()

Solve a system of odes with scipy  how to reference different indexes?
Date : March 29 2020, 07:55 AM
I hope this helps you . The function dfdt takes and returns the state and derivative, respectively as arrays (or other iterables). Thus, all you have to do is to loop over all indices and apply your operations accordingly. For example: def dfdt(t,f):
output = np.empty_like(f)
for i,entry in enumerate(f)
output[i] = f[i]**2 * coeffs[i,:].sum()  2*f[i]*coeffs.sum()
return output
def dfdt(t,f):
return f**2 * coeffs.sum(axis=1)  2*f*coeffs.sum()

Solving a first order system of ODEs using SymPy expressions and SciPy solver
Date : March 29 2020, 07:55 AM
Does that help There is rarely need to use Eq objects in SymPy; equations are best represented by expressions, which are understood to be equated to zero in the context of solvers. So, ode1 and ode2 should just be expressions for the righthand side of the ODE system. Then they can be lambdified together, as follows: ode1 = a*x  b*x*y
ode2 = c*x*y  d*y
F_sympy = sympy.lambdify((x, y, a, b, c, d), [ode1, ode2])
def F_from_sympy(X, t, args):
a, b, c, d = args
x, y = X
return F_sympy(x, y, a, b, c, d)
F_from_sympy([1, 2], np.linspace(0, 1, 100), (3, 4, 5, 6))

Fitting data to system of ODEs using Python via Scipy & Numpy
Date : March 29 2020, 07:55 AM

