  C RUBY-ON-RAILS MYSQL ASP.NET DEVELOPMENT RUBY .NET LINUX SQL-SERVER REGEX WINDOWS ALGORITHM ECLIPSE VISUAL-STUDIO STRING SVN PERFORMANCE APACHE-FLEX UNIT-TESTING SECURITY LINQ UNIX MATH EMAIL OOP LANGUAGE-AGNOSTIC VB6 MSBUILD # Solution from SciPy solve_ivp contains oscillations for a system of first-order ODEs

## Solution from SciPy solve_ivp contains oscillations for a system of first-order ODEs Tag : python , By : Cowtung Date : January 11 2021, 03:34 PM

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

# Calculate Temperatures

interp_qheat = interp1d(time, qheat)

def tc_dt(t, tc, ts, q):
"""
dTc/dt = (Ts-Tc)/(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 = (Tf-Ts)/(Ru*Cs) - (Ts-Tc)/(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 = (Ts-Tc)/(Rc*Cc) + Q/Cc
dTs/dt = (Tf-Ts)/(Ru*Cs) - (Ts-Tc)/(Rc*Cs)

determine the Jacobian matrix of the right-hand side as

Jacobian matrix = [df1/dTc, df2/dTc]
[df1/dTs, df2/dTs]

where

f1 = (Ts-Tc)/(Rc*Cc) + Q/Cc
f2 = (Tf-Ts)/(Ru*Cs) - (Ts-Tc)/(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):
"""
Right-hand side of the system of ODEs.
"""
q = interp_qheat(t)
tcdt = tc_dt(t, y, y, q)
tsdt = ts_dt(t, y, y)
return tcdt, tsdt

t0 = time
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, label='tc')
ax.plot(sol.t, sol.y, 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 .

## Solving a system of odes (with changing constant!) using scipy.integrate.odeint?

Tag : python , By : micaleel
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
y = u
z = u
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()
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 time-dependent 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 time-courses. 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
Zi = y
Ri = y

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')

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()

# 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?

Tag : python , By : ganok_tor
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

Tag : python , By : BooTeK
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 right-hand 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

Tag : python , By : Erik
Date : March 29 2020, 07:55 AM 