Python Forum
How implement parametric ODE equation in python
Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
How implement parametric ODE equation in python
#1
How can I implement following equations in python by using odeint(related to SIS epidemic models). notice that dot up on S means derivation dS/dt and beta, Gamma are constant.

[Image: e4aee24ee052bd1981be4b3b60fd1aac-full.png]

tnx.
Reply
#2
Hi,
I just wrote an example demonstrating the way how to write and solve matrix equations

import numpy as np 
from scipy.integrate import odeint


def equations(y, t, shape=None): # in case of rectangular matrix/equations shape is (n, m), in case of rectangular shape is (n, n); we will pass shape parameter below
    '''
    This defines matrix equations:
      
       dA/dt = A^2 - 4 * A + t,

       where A is a matrix with n rows and n columns 
       (the shape of A is assumed to be n x n). 
    '''
    A = y.reshape(*shape) #lets transform our flat array (y) to a matrix (This allows to avoid elemntwise writing of the equations)
    result = A @ A - 4 * A + t # matrix equations, @ -- denotes matrix multiplication (see PEP 465)
    return result.ravel() # reverse transform (to flat array)


# integration process 

n = 10 # matrix size
A0 = np.random.rand(n, n) # initial conditions, i.e. A0.shape = (n, n)
y0 = A0.ravel() # initial conditions should be flat  (this is needed for odeint, or ode from scipy package)

t = np.linspace(0, 10, 200) # array of points where we are planing to get the solution of the matrix equation.


solution = odeint(equations, y0, t, args=((n, n),)) # args=((n, n),) sets the shape to (n, n) in the equations above

def A(solution, j, shape):
    '''Returns A(t[j])'''
    return solution[j, :].reshape(shape)

print('Result at {} is {}'.format(t[10], A(solution, 10, shape=(n, n))))
However, your equation needs additional clarifications:
Suppose we have matrices S and I. The equation refers to S[s-1, i+1] elements. These elements (probably, boundary conditions) should be defined, e.g. when s=1 this requires S[0, i+1] to be defined, and when i=N (N denotes equation's size) this requires S[s-1, N + 1].

So, you need to define boundary conditions, e.g. S[0, i], S[s, N + 1], I[0, i], I[s, M+1] and shapes of the your matrices (M, N).

if S.shape is NxM and I.shape is NxM, then
1) y.shape will be 2*N*M, y = [elements of S (size=N*M), elements of I (size=N*M)]
2) transformations inside the function that defines equations should be: S = y[:N*M].reshape(N, M), I = y[N*M:].reshape(N, M)
3) you probably will need to rewrite your equations using Cython (scipy.weave, C). If loops in your equations (i.e. double sums) write in pure Python, computational process will be quite slow.
Reply


Possibly Related Threads…
Thread Author Replies Views Last Post
  Python sympy problem for symbolic equation definition joellapointe_engineering 0 87 Mar-24-2024, 11:09 PM
Last Post: joellapointe_engineering
  Parametric portfolio optimization by Brandt 2009 schnellinga 5 2,847 Jul-26-2022, 07:16 PM
Last Post: paulyan

Forum Jump:

User Panel Messages

Announcements
Announcement #1 8/1/2020
Announcement #2 8/2/2020
Announcement #3 8/6/2020