model

model#

model

class pygom.model.CompileCanary[source]#

Hold the need for (re-)compilation for various functions

A subclass of this should specify the states to watch

They may all be tripped to True using the trip() method An individual may be reset with the reset() method or with a direct assignment (they may not be tripped in this way).

reset(name)[source]#

Reset a canary

Parameters: name (string) – the name of the canary to reset
Return type: None

trip()[source]#: Trip all the canaries :rtype: None

pygom.model.DFE(ode, disease_state)[source]#

Returns the disease free equilibrium from an ode object

Parameters

ode (BaseOdeModel) – a class object from pygom
diseaseState (array like) – name of the disease states

Returns

e – disease free equilibrium

Return type

array like

class pygom.model.DeterministicOde(state=None, param=None, derived_param=None, transition=None, event=None, birth_death=None, ode=None)[source]#

This contains the interface and operation built above the already defined set of ode

Parameters

state (list) – A list of states (string)
param (list) – A list of the parameters (string)
derived_param (list) – A list of the derived parameters (tuple of (string,string))
transition (list) – A list of transition (Transition)
birth_death (list) – A list of birth or death process (Transition)
ode (list) – A list of ode (Transition)

add_compiled_sympy_object(method_name, compiled_obj_name, sympy_obj_generator_func, oT, is_master_canary)[source]#: Take sympy_obj_generator_func

add_func(method_name, sympy_obj_generator_func, oT=None, is_master_canary=False)[source]#: Template to add a method, called with method_name, which gives the numerical output of compiled_fn_name, which is the compiled version of the sympy object generated by sympy_obj_generator_func. This carries out the essential check to see if recompilation is necessary and acts accordingly.

adjoint(state, t, state_param, func=None)[source]#

Compute the adjoint given the adjoint vector, time, state variable and the objective function. Note that this function is very restrictive in the sense that the (original) state variable changes through time but this assumes it is a constant, i.e. we assume that the original system is linear.

Parameters

state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.
t (double) – The current time.
state_param (array like) – The state vector that is (or maybe) required to evaluate the jacobian of the original system
func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns

output of the same length as the ode

Return type

numpy.ndarray

Notes

The size of lambda should be the same as the state. The integral should be starting from T, the final time of the original system and is integrated backwards (for stability).

adjoint_T(t, state, state_param, func=None)[source]#: Same as adjoint() but with t as first parameter

adjoint_interpolate(state, t, interpolant, func=None)[source]#

Compute the adjoint given the adjoint vector, time, the functions which was used to interpolate the state variable

Parameters

state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.
t (double) – The current time.
interpolant (list) – list of interpolating functions of the state
func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns

output of the same length as the ode

Return type

numpy.ndarray

adjoint_interpolate_T(t, state, interpolant, objInput=None)[source]#: Same as adjoint_interpolate() but with t as first parameter

adjoint_interpolate_jacobian(state, t, interpolant, func=None)[source]#

Compute the Jacobian of the adjoint given the adjoint vector, time, function of the interpolation on the state variables and the objective function. This is simply the same as the negative Jacobian of the ode transposed.

Parameters

state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.
t (double) – The current time.
interpolant (list) – list of interpolating functions of the state
func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns

output of is a two dimensional array of size [number of state x number of state]

Return type

numpy.ndarray

Notes

Same as adjoint_jacobian() but takes a list of interpolating function instead of a single (vector) value

See also

adjoint_jacobian()

adjoint_interpolate_jacobian_T(t, state, interpolant, func=None)[source]#: Same as adjoint_interpolate_jacobian() but with t as first parameter

adjoint_jacobian(state, t, state_param, func=None)[source]#

Compute the jacobian of the adjoint given the adjoint vector, time, state variable and the objective function. This is simply the same as the negative jacobian of the ode transposed.

Parameters

state (array like) – The current value of lambda, where lambda’s are the Lagrangian multipliers of the differential equation.
t (double) – The current time.
state_param (array like) – The state vector that is (or maybe) required to evaluate the jacobian of the original system
func (callable) – This should take inputs similar to an ode, i.e. of the form func(y,t). If j(y,t) is the cost function, then func is a function that calculates \(\partial j \over \partial x\).

Returns

output of is a two dimensional array of size [number of state x number of state]

Return type

numpy.ndarray

Notes

It takes the same number of argument as the adjoint for simplicity when integrating.

See also

adjoint()

adjoint_jacobian_T(t, state, state_param, func=None)[source]#: Same as adjoint_jacobian_T() but with t being the first parameter

diff_jacobian_T(t, state)[source]#: Same as diff_jacobian() but with t as first parameter

eval_forwardforward(FF, S, state, t)[source]#

Evaluate a single f(x) of the forward-forward sensitivities

Parameters

FF (array like) – this is in fact a 3rd order Tensor, aka 3d array
S (array like) – sensitivities in matrix form
state (array like) – the current state
t (numeric) – time

Returns

f(x) of size [number of state * (number of parameters * number of parameters)]

Return type

numpy.ndarray

eval_hessian(parameters=None, time=None, state=None)[source]#

Evaluate the hessian given parameters, state and time. An extension of hessian() but now also include the parameters.

Parameters

parameters (list) – see parameters()
time (double) – The current time
state (array list) – The current numerical value for the states which can be numpy.ndarray or list

Returns

list of dimension number of state, each with matrix [number of parameters x number of parameters] in sympy.matricies.matricies

Return type

list

See also

grad(), eval_grad()

eval_sens_jacobian_state(time=None, state=None, sens=None)[source]#

Evaluate the jacobian of the sensitivities w.r.t the states given parameters, state and time. An extension of sens_jacobian_state() but now also include the parameters.

Parameters

parameters (list) – see parameters()
time (double) – The current time
state (array list) – The current numerical value for the states which can be numpy.ndarray or list

Returns

Matrix of dimension [number of state x number of state]

Return type

numpy.matrix or mpmath.matrix

Notes

Name and order of state and time are also different.

See also

sens_jacobian_state()

eval_sensitivity(S, t, state, by_state=False)[source]#

Evaluate the sensitivity given state and time

Parameters

S (array like) – Which should be numpy.ndarray. The starting sensitivity of size [number of state x number of parameters]. Which are normally zero or one, depending on whether the initial conditions are also variables.
t (double) – The current time
state (array like) – The current numerical value for the states which can be numpy.ndarray or list
by_state (bool) – how we want the output to be arranged. Default is True so that we have a block diagonal structure

Return type

numpy.ndarray

Notes

It is different to eval_ode() and eval_jacobian() in that the extra input argument is not a parameter

See also

sensitivity()

eval_sensitivityIV(S, IV, t, state)[source]#

Evaluate the sensitivity with initial values given state and time

Parameters

S (array like) – Which should be numpy.ndarray. The starting sensitivity of size [number of state x number of parameters]. Which are normally zero or one, depending on whether the initial conditions are also variables.
IV (array like) – sensitivities for the initial values
t (double) – The current time
state (array like) – The current numerical value for the states which can be numpy.ndarray or list

Returns

\(f(s(x,\theta))\) and \(f(s(x_{0}))\)

Return type

numpy.ndarray

Notes

It is different to eval_ode() and eval_jacobian() in that the extra input argument is not a parameter.

See also

sensitivityIV()

forwardforward(ff, t, state, s)[source]#

Evaluate a single \(f(x)\) of the forward-forward sensitivities

Parameters

ff (array like) – the forward-forward sensitivities in vector form
t (numeric) – time
state (array like) – the current state
s (array like) – forward sensitivities in vector form

Returns

\(f(x)\) of size [number of state * (number of parameters * number of parameters)]

Return type

numpy.ndarray

forwardforward_T(t, ff, s, state)[source]#: Same as forwardforward() but with t as the first parameter

get_diff_jacobian_eqn()[source]#

Returns the jacobian differentiate w.r.t. states in algebraic form

Returns: list of size (num of state,) each with sympy.matrices.matrices of dimension [number of state x number of state]
Return type: list

get_grad_eqn()[source]#

Return the gradient of the ode in algebraic form

Returns: A matrix of dimension [number of state x number of parameters]
Return type: sympy.matrices.matrices

get_grad_jacobian_eqn()[source]#

Return the jacobian of the gradient in algebraic form

Returns: A matrix of dimension [number of state * number of parameters x number of state]
Return type: sympy.matrices.matrices

See also

get_grad_eqn()

get_hessian_eqn()[source]#

Return the Hessian of the ode in algebraic form

Returns: list of dimension number of state, each with matrix [number of parameters x number of parameters] in sympy.matricies.matricies
Return type: list

Notes

We deliberately return a list instead of a 3d array of a tensor to avoid confusion

get_jacobian_eqn()[source]#

Returns the jacobian in algebraic form

Returns: A matrix of dimension [number of state x number of state]
Return type: sympy.matrices.matrices

get_ode_eqn()[source]#: Build the algebraic system of ODE’s given the transitions and events.

get_transition_graph(file_name=None, show=True)[source]#

Returns the transition graph using graphviz

Parameters

file_name (str, optional) – name of the output file, defaults to None
show (bool, optional) – If the graph should be plotted, defaults to True

Return type

graphviz.Digraph

grad_T(t, state)[source]#: Same as grad_T() but with t as first parameter

grad_jacobianT(t, state)[source]#: Same as grad_jacobian() but with t as first parameter

hessian(state, time)[source]#

Evaluate the hessian given state and time

Parameters

state (array like) – The current numerical value for the states which can be numpy.ndarray or list
t (double) – The current time

Returns

list of dimension number of state, each with matrix [number of parameters x number of parameters] in sympy.matricies.matricies

Return type

list

property initial_state#: Return the initial state values

property initial_time#: Return the initial time

property initial_values#: Returns the initial values, both time and state as a tuple (x0, t0)

integrate(t, full_output=False)[source]#

Integrate over a range of t when t is an array and a output at time t

Parameters

t (array like) – the range of time points which we want to see the result of
full_output (bool) – if we want additional information

integrate2(t, full_output=False, method=None)[source]#

Integrate over a range of t when t is an array and a output at time t. Select a suitable method to integrate when method is None.

Parameters

t (array like) – the range of time points which we want to see the result of
full_output (bool) – if we want additional information
method (str, optional) – the integration method. All those available in ode are allowed with ‘vode’ and ‘ivode’ representing the non-stiff and stiff version respectively. Defaults to None, which tries to choose the integration method via eigenvalue analysis (only one) using the initial conditions

is_stiff(state=None, t=None)[source]#

Test on the eigenvalues of the jacobian. We classify the problem as stiff if any of the eigenvalues are positive

Parameters

state (array like) – The current numerical value for the states which can be numpy.ndarray or list
t (double) – The current time

Returns

eigenvalues of the system given input

Return type

numpy.ndarray

jacobian_T(t, state)[source]#: Same as jacobian() but with t as first parameter

jacobian_eigenvalue(state=None, t=None)[source]#

Find out the eigenvalues of the jacobian given state and time. If None is given, the initial values are used.

Parameters

state (array like) – The current numerical value for the states which can be numpy.ndarray or list
t (double) – The current time

Returns

True if any eigenvalue is positive

Return type

bool

linear_ode()[source]#

To check whether the input ode is linear

Returns: True if it is linear, False otherwise
Return type: bool

ode_T(t, state)[source]#: Same as ode() but with t as the first parameter

ode_and_forwardforward(state_param, t)[source]#

Evaluate a single f(x) of the ode and the forward-forward sensitivities

Parameters

state_param (array like) – state and forward-forward sensitivities in vector form
t (numeric) – time

Returns

same size as the state_param input

Return type

numpy.ndarray

ode_and_forwardforward_T(t, state_param)[source]#: Same as odeAndForwardForward() but with time as the first input

ode_and_forwardforward_jacobian(state_param, t)[source]#

Return the jacobian after evaluation given the input of the state and the forward forward sensitivities

Parameters

state_param (array like) – state and forward-forward sensitivities in vector form
t (numeric) – time

Returns

size of (a,a) where a is the length of the state_param input

Return type

numpy.ndarray

ode_and_forwardforward_jacobian_T(t, state_param)[source]#: Same as ode_and_forwardforward_jacobian() but with t being the first parameters

ode_and_sensitivity(state_param, t, by_state=False)[source]#

Evaluate the sensitivity given state and time

Parameters

state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.
t (double) – The current time
by_state (bool) – Whether the output vector should be arranged by state or by parameters. If False, then it means that the vector of output is arranged according to looping i,j from Sensitivity_{i,j} with i being the state and j the param. This is the preferred way because it leds to a block diagonal Jacobian

Returns

concatenation of 2 element. First contains the ode, second the sensitivity. Both are of type numpy.ndarray

Return type

list

See also

sensitivity(), ode()

ode_and_sensitivityIV(state_param, t)[source]#

Evaluate the sensitivity given state and time

Parameters

state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.
t (double) – The current time

Returns

concatenation of 3 element. First contains the ode, second the sensitivity, then the sensitivity of the initial value. All of them are of type numpy.ndarray

Return type

list

See also

sensitivity(), ode()

ode_and_sensitivityIV_T(t, state_param)[source]#: Same as ode_and_sensitivityIV() but with t as first parameter

ode_and_sensitivityIV_jacobian(state_param, t)[source]#

Evaluate the sensitivity given state and time. Output a block diagonal sparse matrix as default.

Parameters

state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.
t (double) – The current time
byState (bool) – How the output is arranged, according to the vector of output. It can be in terms of state or parameters, where by state means that the jacobian is a block diagonal matrix.

Returns

output of a square matrix of size: number of ode + 1 times number of parameters

Return type

numpy.ndarray

See also

ode_and_sensitivity()

ode_and_sensitivityIV_jacobian_T(t, state_param)[source]#: Same as ode_and_sensitivityIV_jacobian() but with t as first parameter

ode_and_sensitivity_T(t, state_param, by_state=False)[source]#: Same as ode_and_sensitivity() but with t as first parameter

ode_and_sensitivity_jacobian(state_param, t, by_state=False)[source]#

Evaluate the sensitivity given state and time. Output a block diagonal sparse matrix as default.

Parameters

state_param (array like) – The current numerical value for the states as well as the sensitivities values all in one. We assume that the state values comes first.
t (double) – The current time
by_state (bool) – How the output is arranged, according to the vector of output. It can be in terms of state or parameters, where by state means that the jacobian is a block diagonal matrix.

Returns

output of a square matrix of size: number of ode + 1 times number of parameters

Return type

numpy.ndarray

See also

ode_and_sensitivity()

ode_and_sensitivity_jacobian_T(t, state_param, by_state=False)[source]#: Same as ode_and_sensitivity_jacobian() but with t as first parameter

plot()[source]#

Plot the results of the integration

Notes

If we have 3 states or more, it will always be arrange such that it has 3 columns. Uses the operation from odeutils

print_ode(latex_output=False)[source]#

Prints the ode in symbolic form onto the screen/console in actual symbols rather than the word of the symbol.

Parameters: latex_output (bool, optional) – Defaults to false which prints the equation in terms of symbols, if set to yes then the formula in terms of latex equations will be printed onto the screen.

sens_jacobian_state(state_param, t)[source]#

Evaluate the jacobian of the sensitivity w.r.t. the state given state and time

Parameters

state_param (array like) – The current numerical value for the states as well as the sensitivities, which can be numpy.ndarray or list
t (double) – The current time

Returns

Matrix of dimension [number of state * number of parameters x number of state]

Return type

numpy.ndarray

sens_jacobian_state_T(t, state)[source]#: Same as sens_jacobian_state_T() but with t as first parameter

sensitivity(sens, t, state, by_state=False)[source]#

Evaluate the sensitivity given state and time. The default is to output the values by parameters, i.e. \(s_{i},\ldots,s_{i+n}\) are partial derivatives w.r.t. the states for \(i \in {1,1+p,1+2p,1+3p, \ldots, 1+(n-1)p}\). This is to take advantage of the fact that we have a block diagonal jacobian that was already evaluated

Parameters

sens (array like) – The starting sensitivity of size [number of state x number of parameters]. Which are normally zero or one, depending on whether the initial conditions are also variables.
t (double) – The current time
state (array like) – The current numerical value for the states which can be numpy.ndarray or list
by_state (bool) – how we want the output to be arranged. Default is True so that we have a block diagonal structure

Return type

numpy.ndarray

sensitivityIV(sensIV, t, state)[source]#

Evaluate the sensitivity which include the initial values as our parameters given state and time. The default is to output the values by parameters, i.e. \(s_{i},\ldots,s_{i+n}\) are partial derivatives w.r.t. the states for \(i \in {1,1+p,1+2p,1+3p, \ldots, 1+(n-1)p}\). This is to take advantage of the fact that we have a block diagonal Jacobian that was already evaluated.

Parameters

sensIV (array like) – The starting sensitivity of size [number of state x number of parameters] + [number of state x number of state] for the initial condition. The latter is an identity matrix at time zero.
t (double) – The current time
state (array like) – The current numerical value for the states which can be numpy.ndarray or list

Returns

output of the same length as the ode

Return type

numpy.ndarray

sensitivityIV_T(t, sensIV, state)[source]#: Same as sensitivityIV() but with t as first parameter

sensitivity_T(t, sens, state, by_state=False)[source]#: Same as sensitivity() but with t as first parameter

class pygom.model.Event(transition_list, rate=None)[source]#: Class to contain transitions

pygom.model.FitzHugh(param=None)[source]#

The standard FitzHugh model without external input [FitzHugh1961]

\[\begin{split}\frac{dV}{dt} &= c ( V - \frac{V^{3}}{3} + R) \\ \frac{dR}{dt} &= -\frac{1}{c}(V - a + bR).\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.FitzHugh({'a':0.2, 'b':0.2, 'c':3.0})
>>> t = np.linspace(0, 20, 101)
>>> x0 = [1.0, -1.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

pygom.model.Influenza_SLIARD(param=None)[source]#: A simple influenza model from [Brauer2008], page 323.

\[\begin{split}\frac{dS}{dt} &= -S \beta (I + \delta A) \\ \frac{dL}{dt} &= S \beta (I + \delta A) - \kappa L \\ \frac{dI}{dt} &= p \kappa L - \alpha I \\ \frac{dA}{dt} &= (1 - p) \kappa L - \eta A \\ \frac{dR}{dt} &= f \alpha I + \eta A \\ \frac{dD}{dt} &= (1 - f) \alpha I\end{split}\]

pygom.model.Legrand_Ebola_SEIHFR(param=None)[source]#

The Legrand Ebola model [Legrand2007] with 6 compartments that includes the H = hospitalization and F = funeral state. Note that because this is an non-autonomous system, there are in fact a total of 7 states after conversion. The set of equations that describes the model are

\[\begin{split}\frac{dS}{dt} &= -(\beta_{I}SI + \beta_{H}SH + \beta_{F}SF) \\ \frac{dE}{dt} &= (\beta_{I}SI + \beta_{H}SH + \beta_{F}SF) - \alpha E \\ \frac{dI}{dt} &= \alpha E - (\gamma_{H} \theta_{1} + \gamma_{I}(1-\theta_{1})(1-\delta_{1}) + \gamma_{D}(1-\theta_{1})\delta_{1})I \\ \frac{dH}{dt} &= \gamma_{H}\theta_{1}I - (\gamma_{DH}\delta_{2} + \gamma_{IH}(1-\delta_{2}))H \\ \frac{dF}{dt} &= \gamma_{D}(1-\theta_{1})\delta_{1}I + \gamma_{DH}\delta_{2}H - \gamma_{F}F \\ \frac{dR}{dt} &= \gamma_{I}(1-\theta_{1})(1-\delta_{1})I + \gamma_{IH}(1-\delta_{2})H + \gamma_{F}F.\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> x0 = [1.0, 3.0/200000.0, 0.0, 0.0, 0.0, 0.0, 0.0]
>>> t = np.linspace(0, 25, 100)
>>> ode = common_models.Legrand_Ebola_SEIHFR([('beta_I',0.588),('beta_H',0.794),('beta_F',7.653),('omega_I',10.0/7.0),('omega_D',9.6/7.0),('omega_H',5.0/7.0),('omega_F',2.0/7.0),('alphaInv',7.0/7.0),('delta',0.81),('theta',0.80),('kappa',300.0),('interventionTime',7.0)])
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

pygom.model.Lorenz(param=None)[source]#

Lorenz attractor define by three parameters, \(\beta,\sigma,\rho\) as per [Lorenz1963].

\[\begin{split}\frac{dx}{dt} &= \sigma (y-x) \\ \frac{dy}{dt} &= x (\rho - z) - y \\ \frac{dz}{dt} &= xy - \beta z\end{split}\]

Examples

>>> import matplotlib.pyplot as plt
>>> import numpy
>>> from pygom import common_models
>>> t = numpy.linspace(0, 20, 101)
>>> params = {'beta':8.0/3.0, 'sigma':10.0, 'rho':28.0}
>>> ode = common_models.Lorenz(params)
>>> ode.initial_values = ([1., 1., 1.], t[0])
>>> solution = ode.integrate(t[1::])
>>> plt.plot(solution[:,0], solution[:,2])
>>> plt.show()

pygom.model.Lotka_Volterra(param=None)[source]#

Standard Lotka-Volterra model with two states and four parameters [Lotka1920]

\[\begin{split}\frac{dx}{dt} &= \alpha x - cxy \\ \frac{dy}{dt} &= -\delta y + \gamma xy\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> params = {'alpha':1, 'delta':3, 'c':2, 'gamma':6}
>>> ode = common_models.Lotka_Volterra(params)
>>> ode.initial_values = ([2.0, 6.0], 0)
>>> t = np.linspace(0.1, 100, 10000)
>>> ode.integrate(t)
>>> ode.plot()

class pygom.model.ODEVariable(ID, name=None, units=None, real=True)[source]#

A class that defines the variables in our ODE

Parameters

ID (str) – identifier of the variable
name (str, optional) – name of the variable in human readable format. Defaults to None, which then takes the ID as the name
units (str, optional) – what unit the variable takes. Defaults to None.
real (bool, optional) – if the variable can only be a real number, defaults to True

class pygom.model.OrderedDict[source]#

Dictionary that remembers insertion order

clear() → None. Remove all items from od.#

copy() → a shallow copy of od#

fromkeys(value=None)#: Create a new ordered dictionary with keys from iterable and values set to value.

items() → a set-like object providing a view on D's items#

keys() → a set-like object providing a view on D's keys#

move_to_end(key, last=True)#

Move an existing element to the end (or beginning if last is false).

Raise KeyError if the element does not exist.

pop(key[, default]) → v, remove specified key and return the corresponding value.#: If the key is not found, return the default if given; otherwise, raise a KeyError.

popitem(last=True)#

Remove and return a (key, value) pair from the dictionary.

Pairs are returned in LIFO order if last is true or FIFO order if false.

setdefault(key, default=None)#

Insert key with a value of default if key is not in the dictionary.

Return the value for key if key is in the dictionary, else default.

update([E, ]**F) → None. Update D from dict/iterable E and F.#: If E is present and has a .keys() method, then does: for k in E: D[k] = E[k] If E is present and lacks a .keys() method, then does: for k, v in E: D[k] = v In either case, this is followed by: for k in F: D[k] = F[k]

values() → an object providing a view on D's values#

pygom.model.R0(ode, disease_state)[source]#

Returns the basic reproduction number, in symbolic form when the parameter values are not available

Parameters

ode (BaseOdeModel) – a class object from pygom
diseaseStateIndex (array like) – name of the disease states

Returns

e – R0

Return type

array like

See also

getDiseaseProgressionMatrices(), getR0GivenMatrix()

pygom.model.R0_from_matrix(F, V, disease_state=None)[source]#

Returns the symbolic form of the basic reproduction number. This will include the states symbols which is different from getR0() where the states is replaced by the values of the disease-free equilibrium.

Parameters

F (sympy.matrices.MatrixBase) – secondary infection rates
V (sympy.matrices.MatrixBase) – disease progression rates
disease_state (list like, optional) – list of the disease state as sympy.Symbol. Defaults to None which assumes that \(F,V\) had been differentiated

Returns

e – the eigenvalues of \(FV^{-1}\) for the disease states

Return type

sympy.matrices.MatrixBase

See also

getDiseaseProgressionMatrices(), getR0()

pygom.model.Robertson(param=None)[source]#

The so called Robertson problem [Robertson1966], which is a standard example used to test stiff integrator.

\[\begin{split}\frac{dy_{1}}{dt} &= -0.04 y_{1} + 1 \cdot 10^{4} y_{2} y_{3} \\ \frac{dy_{2}}{dt} &= 0.04 y_{1} - 1 \cdot 10^{4} y_{2} y_{3} - 3 \cdot 10^{7} y_{2}^{2}\\ \frac{dy_{3}}{dt} &= 3 \cdot 10^{7} y_{2}^{2}\end{split}\]

Examples

>>> from pygom import common_models
>>> import numpy
>>> t = numpy.append(0, 4*numpy.logspace(-6, 6, 1000))
>>> ode = common_models.Robertson()
>>> ode.initial_values = ([1.0,0.0,0.0], t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot() # note that this is not being plotted in the log scale

pygom.model.SEIR(param=None, init=None)[source]#: A standard SIR model [Brauer2008] with population N

pygom.model.SEIR_Birth_Death(param=None)[source]#

A standard SEIR model with birth and death [Aron1984], defined by the ode

\[\begin{split}\frac{dS}{dt} &= \mu - \beta SI - \mu S \\ \frac{dE}{dt} &= \beta SI - (\mu + \alpha) E \\ \frac{dI}{dt} &= \alpha E - (\mu + \gamma) I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

Uses the same set of parameters as the examples in SEIR() apart from \(\mu\) which is new.

>>> import numpy as np
>>> from pygom import common_models
>>> params = {'beta':1800, 'gamma':100, 'alpha':35.84, 'mu':0.02}
>>> ode = common_models.SEIR_Birth_Death(params)
>>> t = np.linspace(0, 50, 1001)
>>> x0 = [0.0658, 0.0007, 0.0002, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution,output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()

See also

SEIR()

pygom.model.SEIR_Birth_Death_Periodic(param=None)[source]#

A SEIR birth death model with periodic contact [Aron1984], defined by the ode

\[\begin{split}\frac{dS}{dt} &= \mu - \beta(t)SI - \mu S \\ \frac{dE}{dt} &= \beta(t)SI - (\mu + \alpha) E \\ \frac{dI}{dt} &= \alpha E - (\mu + \gamma) I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

where

\[\beta(t) = \beta_{0} (1 + \beta_{1} \cos(2 \pi t)).\]

An extension of an SEIR birth death model by varying the contact rate through time.

Examples

Uses the same set of parameters as the examples in SEIR_Birth_Death() but now we have two beta parameters instead of one.

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from pygom import common_models
>>> params = {'beta0':1800, 'beta1':0.2, 'gamma':100, 'alpha':35.84, 'mu':0.02}
>>> ode = common_models.SEIR_Birth_Death_Periodic(params)
>>> t = np.linspace(0, 50, 1001)
>>> x0 = [0.0658, 0.0007, 0.0002, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution,output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()
>>> plt.plot(np.log(solution[:,0]), np.log(solution[:,1]))
>>> plt.show()
>>> plt.plot(np.log(solution[:,0]), np.log(solution[:,2]))
>>> plt.show()

See also

SEIR(), SEIR_Birth_Death(), SIR_Periodic()

pygom.model.SEIR_Birth_Death_Periodic_Waning_Intro(param=None)[source]#: SEIR model with vital dynamics, periodic infectivity, immune waning and external introductions

pygom.model.SEIR_Multiple(n=2, param=None)[source]#

An SEIR model that describe spatial heterogeneity [Brauer2008], page 180. The model originated from [Lloyd1996] and notations used here follows [Brauer2008].

\[\begin{split}\frac{dS_{i}}{dt} &= dN_{i} - dS_{i} - \lambda_{i} S_{i} \\ \frac{dE_{i}}{dt} &= \lambda_{i}S_{i} - (d + \epsilon) E_{i} \\ \frac{dI_{i}}{dt} &= \epsilon E_{i} - (d + \gamma) I_{i} \\ \frac{dR_{i}}{dt} &= \gamma I_{i} - dR_{i}\end{split}\]

where

\[\lambda_{i} = \sum_{j=1}^{n} \beta_{i,j} I_{j} (1\{i \neq j\} p)\]

with \(n\) being the number of patch and \(p\) the coupled factor.

Examples

Use the initial conditions that were derived from the stationary condition specified in [Brauer2008].

>>> import numpy as np
>>> from pygom import common_models
>>> paramEval = {'beta_00':0.0010107, 'beta_01':0.0010107,
>>>              'beta_10':0.0010107, 'beta_11':0.0010107,
>>>              'd':0.02,'epsilon':45.6, 'gamma':73.0,
>>>              'N_0':10**6,'N_1':10**6,'p':0.01}
>>> x0 = [36139.3224081278, 422.560577637822,
>>>       263.883351688369, 963174.233662546]
>>> ode = common_models.SEIR_Multiple()
>>> t = np.linspace(0, 40, 100)
>>> x01 = []
>>> for s in x0:
>>>     x01 += [s]
>>>     x01 += [s]
>>> ode.parameters = paramEval
>>> ode.initial_values = (x01, t[0])
>>> solution, output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()

pygom.model.SIR(param=None)[source]#

A standard SIR model as per [Brauer2008]

\[\begin{split}\frac{dS}{dt} &= -\ \frac{beta SI}{N} \\ \frac{dI}{dt} &= \ \frac{beta SI}{N} - \gamma I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

The model that produced top two graph in Figure 1.3 of the reference above. First, when everyone is susceptible and only one individual was infected.

>>> import numpy as np
>>> from pygom import common_models
>>> N=1e5
>>> ode = common_models.SIR({'beta':0.5, 'gamma':0.2, 'N':N})
>>> t = np.linspace(0, 730, 1001)
>>> i0=1
>>> x0 = [N-i0, i0, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

pygom.model.SIR_Birth_Death(param=None)[source]#

Extension of the standard SIR model [Brauer2008] to also include birth and death

\[\begin{split}\frac{dS}{dt} &= B -\beta SI - \mu S \\ \frac{dI}{dt} &= \beta SI - \gamma I - \mu I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

The model that produced bottom graph in Figure 1.3 of the reference above.

>>> import numpy as np
>>> from pygom import common_models
>>> B = 126372.0/365.0
>>> N = 7781984.0
>>> params = {'beta':3.6, 'gamma':0.2, 'B':B/N, 'mu':B/N}
>>> ode = common_models.SIR_Birth_Death(params)
>>> t = np.linspace(0, 35*365, 10001)
>>> x0 = [0.065, 123.0*(5.0/30.0)/N, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution,output = ode.integrate(t[1::], full_output=True)
>>> ode.plot()

See also

SIR()

pygom.model.SIR_norm(param=None)[source]#

A normalized SIR model:

\[\begin{split}\frac{dS}{dt} &= -\beta SI \\ \frac{dI}{dt} &= \beta SI - \gamma I \\ \frac{dR}{dt} &= \gamma I\end{split}\]

Examples

The model that produced top two graph in Figure 1.3 of the reference above. First, when everyone is susceptible and only one individual was infected.

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.SIR({'beta':3.6, 'gamma':0.2})
>>> t = np.linspace(0, 730, 1001)
>>> N = 7781984.0
>>> x0 = [1.0, 10/N, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

Second model with a more realistic scenario

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.SIR({'beta':3.6, 'gamma':0.2})
>>> t = np.linspace(0, 730, 1001)
>>> N = 7781984.0
>>> x0 = [0.065, 123*(5.0/30.0)/N, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

pygom.model.SIS(param=None)[source]#

Susceptible Infected Susceptible model

\[\begin{split}\frac{dS}{dt} &= -\beta SI + \gamma I \\ \frac{dI}{dt} &= \beta SI - \gamma I\end{split}\]

Examples

>>> import numpy as np
>>> from pygom import common_models
>>> ode = common_models.SIS({'beta':0.5, 'gamma':0.2})
>>> t = np.linspace(0, 20, 100)
>>> x0 = [1.0, 0.1]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

pygom.model.SIS_Periodic(param=None)[source]#

A SIS model with periodic contact, defined by the ode as per [Hethcote1973]

\[\frac{dI}{dt} = (\beta(t)N - \alpha) I - \beta(t)I^{2}\]

where

\[\beta(t) = 2 - 1.8 \cos(5t).\]

As the name suggests, it achieves a (stable) periodic solution.

Examples

>>> from pygom import common_models
>>> import numpy as np
>>> ode = common_models.SIS_Periodic({'alpha':1.0})
>>> t = np.linspace(0, 10, 101)
>>> x0 = [0.1, 0.0]
>>> ode.initial_values = (x0, t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

class pygom.model.SimulateOde(state=None, param=None, derived_param=None, transition=None, event=None, birth_death=None, ode=None)[source]#

This builds on top of DeterministicOde which we simulate the outcome instead of solving it deterministically

Parameters

state (list) – A list of states (string) or (string, (numeric, numeric)) if specifying limits
param (list) – A list of the parameters (string)
derived_param (list) – A list of the derived parameters (tuple of (string,string))
transition (list) – A list of transition (Transition) #TODO Now this might actually only be deterministic ODE objects. Check.
transition – A list of events (Event)
birth_death (list #TODO Now these are wrapped in events. Should try to make work for back compat.) – A list of birth or death process (Transition)
ode (list) – A list of ode (Transition)

cle(x0, t0, t1, output_time=False)[source]#

Stochastic simulation using the CLE approximation starting from time t0 to t1 with the starting state values of x0. The CLE approximation is performed using a simple Euler-Maruyama method with step size h. We assume that the input parameter transition_func provides \(f(x,t)\) while the CLE is defined as \(dx = x + V*h*f(x,t) + \sqrt(f(x,t))*Z*\sqrt(h)\) with \(Z\) being standard normal random variables.

Parameters

x (array like) – state vector
t0 (double) – start time
t1 (double) – final time

exact(x0, t0, t1, output_time=False)[source]#

Stochastic simulation using an exact method starting from time t0 to t1 with the starting state values of x0

Parameters

x (array like) – state vector
t0 (double) – start time
t1 (double) – final time

get_TransitionJacobian()[source]#: Evaluate equation (7) from https://people.cs.vt.edu/~ycao/publication/newstepsize.pdf where F_[i,j] is the change in transition rate a[i], if a transition of type j occurs: F_[i,j] = sum_k diff(a[i], x_k) v_[k,j] where k=state and v[k,j] is how much state x_k changes by if transition of type j occurs.

get_TransitionMean()[source]#: This is the mean and variance of the changes in the transition rates (aka propensity funtions) after a potential timestep: equations (8a) and (8b) from https://people.cs.vt.edu/~ycao/publication/newstepsize.pdf For n transitions the outputs are 2 vectors, each of length n. Outputs are added to self as mu and sigma2

get_TransitionVar()[source]#: This is the mean and variance of the changes in the transition rates (aka propensity funtions) after a potential timestep: equations (8a) and (8b) from https://people.cs.vt.edu/~ycao/publication/newstepsize.pdf For n transitions the outputs are 2 vectors, each of length n. Outputs are added to self as mu and sigma2

get_bd_from_ode(A=None)[source]#: Returns a list of:class:Transition from this object by unrolling the odes. All the elements are of TransitionType.B or TransitionType.D

get_transitions_from_ode()[source]#: Returns a list of Transition from this object by unrolling the odes. All the elements are of TransitionType.T

get_unrolled_obj()[source]#: Returns a SimulateOde with the same state and parameters as the current object but with the equations defined by a set of transitions and birth death process instead of say, odes

hybrid(x0, t0, t1, output_time=False)[source]#

Stochastic simulation using an hybrid method that uses either the first reaction method or the \(\tau\)-leap depending on the size of the states and transition rates. Starting from time t0 to t1 with the starting state values of x0.

Parameters

x (array like) – state vector
t0 (double) – start time
t1 (double) – final time

plot(sim_X=None, sim_T=None)[source]#

Plot the results of a simulation

Takes the output of a function like solve_stochast

Parameters

sim_X (list) – of length iteration each with (len(t),len(state)) if t is a vector, else it outputs unequal shape that was record of all the jumps
sim_T (list or numpy.ndarray) – if t is a single value, it outputs unequal shape that was record of all the jumps. if t is a vector, it outputs t so that it is a numpy.ndarray instead

Notes

If either sim_X or sim_T are None the this function will attempt to plot the deterministic ODE

If we have 3 states or more, it will always be arrange such that it has 3 columns. Uses the operation from odeutils

simulate_param(t, iteration, parallel=False, full_output=False)[source]#

Simulate the ode by generating new realization of the stochastic parameters and integrate the system deterministically.

Parameters

t (array like) – the range of time points which we want to see the result of
iteration (int) – number of iterations you wish to simulate
parallel (bool, optional) – Defaults to True
full_output (bool, optional) – if we want additional information, Y_all in the return, defaults to false

Returns

Y (numpy.ndarray) – of shape (len(t), len(state)), mean of all the simulation
Y_all (np.ndarray) – of shape (iteration, len(t), len(state))

solve_determ(t, iteration=None, parallel=False, full_output=False)[source]#

Simulate the ode by generating new realization of the stochastic parameters and integrate the system deterministically.

Parameters

t (array like) – the range of time points which we want to see the result of
iteration (int) – number of iterations you wish to simulate
parallel (bool, optional) – Defaults to True
full_output (bool, optional) – if we want additional information, Y_all in the return, defaults to false

Returns

Y (numpy.ndarray) – of shape (len(t), len(state)), mean of all the simulation
Y_all (np.ndarray) – of shape (iteration, len(t), len(state))

solve_stochast(t, iteration, parallel=False, exact=False, full_output=False)[source]#

Simulate the ode using stochastic simulation. It switches between a first reaction method and a \(\tau\)-leap algorithm internally. When a parallel backend exists, then a new random state (seed) will be used for each processor. This is due to a lack of appropriate parallel seed random number generator in python.

Parameters

t (array like) – the range of time points which we want to see the result of or the final time point
iteration (int) – number of iterations you wish to simulate
parallel (bool, optional) – Defaults to True
exact (bool, optional) – True if exact simulation is desired, defaults to False
full_output (bool, optional) – if we want additional information, sim_T

Returns

sim_X (list) – of length iteration each with (len(t),len(state)) if t is a vector, else it outputs unequal shape that was record of all the jumps
sim_Jump (list of numpy.ndarray) – Number times each transition happens per timestep
sim_T (list or numpy.ndarray) – if t is a single value, it outputs unequal shape that was record of all the jumps. if t is a vector, it outputs t so that it is a numpy.ndarray instead

total_transition(state, t)[source]#

Evaluate the total transition rate given state and time

Parameters

state (array like) – The current numerical value for the states which can be numpy.ndarray or list
t (double) – The current time

Returns

total rate

Return type

float

class pygom.model.Transition(origin=None, equation=None, transition_type='ODE', destination=None, magnitude='1', ID=None, name=None)[source]#

This class carries the information for transitions defined for an ode, which includes the ode itself, a birth death process where only one state is involved and also a transition between two states

Parameters

origin (str) – Origin state.
equation (str) – Equation defining the transition
transition_type (enum or str, optional) – of type TransitionType or one of (‘ODE’, ‘T’, ‘B’, ‘D’) defaults to ‘ODE’
destination (str, optional) – Destination State. If the transition is not between state, such as a birth or death process, then this is is not required. If it is stated as a birth, death or an ode then it throws an error

property destination#

Return the destination state

Returns: The destination state
Return type: string

property equation#

Return the transition _equation

Returns: The transition _equation
Return type: string

is_between_state()[source]#

Return whether it is a transition between two state

Returns: True if it is a transition between two state False if it is only related to the origin state
Return type: bool

property origin#

Return the origin state

Returns: The origin state
Return type: string

property stochastic#

Return the secondary effects

Returns: The destination state
Return type: string

property transition_type#

Return the type of transition

Returns: One of the four type available from transition_type
Return type: transition_type

class pygom.model.TransitionType(value)[source]#

This is an Enum describing the four feasible type of transitions use to define the ode model BaseOdeModel

The following four types of transitions are available.

B = Birth process

D = Death process

T = Transition between states

ODE = ODE _equation

pygom.model.check_array_type(x, accept_booleans=False)[source]#

Check to see if the type of input is suitable. Only operate on one or two dimension arrays

Parameters

x (array like) – which can be either a numpy.ndarray or list or tuple
accept_boolean (boolean) – If true boolean elements are accepted, else they are not.

Returns

x – checked and converted array

Return type

numpy.ndarray

pygom.model.check_dimension(x, y)[source]#

Compare the length of two array like objects. Converting both to a numpy array in the process if they are not already one.

Parameters

x (array like) – first array
y (array like) – second array

Returns

x (numpy.array) – checked and converted first array
y (numpy.array) – checked and converted second array

class pygom.model.compileCode(backend=None)[source]#

A class that compiles an algebraic expression in sympy to a faster numerical file using the appropriate backend.

compileExpr(inputSymb, inputExpr, backend=None, compileType=False)[source]#

Compiles the expression given the symbols. Determines the backend if required.

Parameters

inputSymb (list) – the set of symbols for the input expression
inputExpr (expr) – expression in sympy
backend (optional) – the backend we want to use to compile
compileType (optional) – defaults to False. If True, return an extra output that informs the end user of the method used to compile the equation, can be one of (np, mpmath, sympy)

Return type

Compiled function taking arguments of the input symbols

compileExprAndFormat(inputSymb, inputExpr, backend=None, modules=None, outType=None)[source]#

Compiles the expression given the symbols and determine which type of output is it. Transforms the output appropriately into numpy

Parameters

inputSymb (list) – the set of symbols for the input expression
inputExpr (expr) – expression in sympy
backend (optional) – the backend we want to use to compile
modules (optional) – in the event that f2py and Cython fails, which modules do we want to try and compile against

Return type

Function determined from the input using closures.

pygom.model.disease_progression_matrices(ode, disease_state, diff=True)[source]#

Returns (F,V), the secondary infection rates and disease progression rate respectively.

Parameters

ode (BaseOdeModel) – an ode class in pygom
diseaseStates (array like) – the name of the disease states
diff (bool, optional) – if the first derivative of the matrices are return, defaults to true

Returns

(F, V) – The progression matrices. If diff=False, then we return the \(F_{i}\) and \(V_{i}\) matrices as per [Brauer2008].

Return type

tuple

pygom.model.integrate(ode, x0, t, full_output=False)[source]#

A wrapper on top of odeint using DeterministicOde.

Parameters

ode (object) – of type DeterministicOde
t (array like) – the time points including initial time
full_output (bool, optional) – If the additional information from the integration is required

pygom.model.integrateFuncJac(func, jac, x0, t0, t, args=(), includeOrigin=False, full_output=False, method=None, nsteps=10000)[source]#

A replacement for scipy.integrate.odeint which performs integration using scipy.integrate.ode, tries to pick the correct integration method at the start through eigenvalue analysis

Parameters

func (callable) – the ode \(f(x)\)
jac (callable) – jacobian of the ode, \(J_{i,j} = \nabla_{x_{j}} f_{i}(x)\)
x0 (numpy.ndarray or list of numeric) – initial value of the states
t0 (float) – initial time
args (tuple, optional) – additional arguments to be passed on
includeOrigin (bool, optional) – if the output should include the initial states x0
full_output (bool, optional) – if additional output is required
method (str, optional) – the integration method. All those availble in ode are allowed with ‘vode’ and ‘ivode’ representing the non-stiff and stiff version respectively. Defaults to None, which tries to choose the integration method via eigenvalue analysis (only one) using the initial conditions
nstep (int, optional) – number of steps allowed between each time point of the integration

Returns

solution (array like) – a np.ndarray of shape (len(t), len(x0)) if includeOrigin is False, else an extra row with x0 being the first.

output (dict, only returned if full_output=True) – Dictionary containing additional output information

key	meaning
’ev’	vector of eigenvalues at each t
’maxev’	maximum eigenvalue at each t
’minev’	minimum eigenvalue at each t
’suc’	list whether integration is successful
’in’	name of integrator

pygom.model.is_list_like(x)[source]#

Test whether the input is a type that behaves like a list, such as (list,tuple,np.ndarray)

Parameters: x – anything
Returns: True if it belongs to one of the three expected type (list,tuple,np.ndarray)
Return type: bool

pygom.model.plot_det(solution, t, stateList=None, y=None, yStateList=None)[source]#

Plot the results of the integration

Parameters

solution (numpy.ndarray) – solution from the integration
t (array like) – the vector of time where the integration output correspond to
stateList (list) – name of the states, if available
Notes –
more (If we have 5 states or) –
such (it will always be arrange) –
columns. (that it has 3) –

pygom.model.plot_stoc(solution, t, stochastic_model)[source]#

Plot the results of a stocastic simulation

Parameters

solution – results of the stochastic simulation
t (array like) – the vector of time where the integration output correspond to
stochastic_model – the model from which this simulation was generated
Notes –
more (If we have 5 states or) –
such (it will always be arrange) –
columns. (that it has 3) –

class pygom.model.shapeAdjust(numState, numParam, numTarget=None)[source]#

A class that change vector into matrix and vice versa for vectors used in DeterministicOde

Parameters

numState (int) – number of states
numParam (int) – number of parameters
numTarget (int, optional) – number of targeted states, default assumes that this is the same as numState

kronParam(A, pre=False)[source]#

A sparse multiplication with an identity matrix of size equal to the number of parameters as initialized

Parameters

A (array like) – a 2d array
pre (bool, optional) – If True, then returns \(I \otimes A\). If False then \(A \otimes I\), where \(A\) is the input matrix, \(I\) is the identity matrix and \(\otimes\) is the kron operator

kronState(A, pre=False)[source]#

A sparse multiplication with an identity matrix of size equal to the number of state as initialized

Parameters

A (array like) – a 2d array
pre (bool, optional) – If True, then returns \(I \otimes A\). If False then \(A \otimes I\), where \(A\) is the input matrix, \(I\) is the identity matrix and \(\otimes\) is the kron operator

matToVecFF(FF)[source]#: Transforms the forward forward sensitivity matrix to a vector

matToVecSens(S)[source]#: Transforms the sensitivity matrix to a vector

vecToMatFF(ff)[source]#: Transforms the forward forward sensitivity vector to a matrix

vecToMatSens(s)[source]#: Transforms the sensitivity vector to a matrix

pygom.model.str_or_list(x)[source]#

Test to see whether input is a string or a list. If it is a string, then we convert it to a list.

Parameters: x – str or list
Returns: x in list form
Return type: x

pygom.model.vanDerPol(param=None)[source]#

The van der Pol equation [vanderpol1926], a second order ode

\[y^{\prime\prime} - \mu (1-y^{2}) y^{\prime} + y = 0\]

where \(\mu > 0\). This can be converted to a first order ode by equating \(x = y^{\prime}\)

\[x^{\prime} - \mu (1 - y^{2}) x + y = 0\]

which result in a coupled ode

\[\begin{split}x^{\prime} &= \mu (1 - y^{2}) x - y \\ y^{\prime} &= x\end{split}\]

and this can be solved via standard method.

Examples

>>> from pygom import common_models
>>> import numpy
>>> t = numpy.linspace(0, 20, 1000)
>>> ode = common_models.vanDerPol({'mu':1.0})
>>> ode.initial_values = ([2.0,0.0], t[0])
>>> solution = ode.integrate(t[1::])
>>> ode.plot()

model

Contents

model#