CODEWITHSUNDEEP
scipysympyode
kemal
kemal

Reputation: 33

How to use a sympy generated Jacobi matrix in the solution of an ODE system ?

I have a first order ODE system which is composed of 3 diff. eqn's. I want to to solve it with scipy.integrate.solve_ivp's BDF method. So I need to calculate jacobi matrix of system (and made it with the help of SymPy).

If i didn't misunderstand; according to the scipy.integrate.solve_ivp document, you must introduce jacobien matrix in the form of jac(t,u) where u should be the state variables of your ODE system. To this end i lambdify jacobien matrix properly.

And my problem arises here. Although I am able to calculate jac(t,u) with some (t,u) such as ((1/800),(150,1E-6,3)), I couldn't send array arguments to my jac. when i introduce jac(t,u) as an argument to solve_ivp it gives an error message. So how should i introduce jac matrix ? Or is my lambdify not proper ?

This is my code. Any help i appreciate it. 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def cvs(t,u):
    u1,u2,u3 = u
    def Qmi(t):
        return t**2
    u1p = Qmi(t)*u3
    u2p = (u1**2)*np.cos(2*np.pi*200*t)
    u3p = (np.sin(2*np.pi*t))*u2**-1
    return [u1p,u2p,u3p]
def jac_func():
    ######### DEFINE THE ODE SYSTEM #########
    import sympy
    sympy.init_printing()
    t = sympy.symbols("t")
    Q_mi = sympy.Function("Q_mi")(t)
    u1 = sympy.Function("u1")(t)
    u2 = sympy.Function("u2")(t)
    u3 = sympy.Function("u3")(t)
    Q_mi = t**2
    u1p = (u3*Q_mi)
    u2p = (u1**2)*sympy.cos(2*sympy.pi*200*t)
    u3p = sympy.sin(2*sympy.pi*5*t)*u2**-1
    ####### CALCULATE THE JACOBIEN ########
    ode_rhs = sympy.Matrix([u1p,u2p,u3p])
    ode_var = sympy.Matrix([u1,u2,u3])
    jac = sympy.Matrix([[ode.diff(var) for var in ode_var]for ode in ode_rhs])
    u = (u1,u2,u3)
    jac_np = sympy.lambdify((t,u),jac,"numpy")
    return jac_np

jac_np = jac_func()
U_0 = [500,20,20]
t = np.linspace(0,100,10000)

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

error message:

ValueError                                Traceback (most recent call last)
<ipython-input-1-8b86ffb3a7cf> in <module>()
41 t = np.linspace(0,100,10000)
42 
---> 43 solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

<lambdifygenerated-1> in _lambdifygenerated(t, _Dummy_188)
  1 def _lambdifygenerated(t, _Dummy_188):
  2     [_Dummy_185, _Dummy_186, _Dummy_187] = _Dummy_188
----> 3     return (array([[0, 0, t**2], [2*_Dummy_185*cos(400*pi*t), 0, 0], [0, -sin(10*pi*t)/_Dummy_186**2, 0]]))

ValueError: setting an array element with a sequence.

Upvotes: 0

Views: 966

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

You are getting the problem because you do what the error message says, you are passing an array where the procedure expects a single number. In 

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

you are trying the constant matrix jac_np(t,U_0) to the Jacobian argument. However, at that point t contains all the t values that you want output samples from. A list of [ array, scalar, scalar ] is not compatible with the numpy arrays.

Long story short, remove the arguments, pass the Jacobian as callable function, as you quite probably intended,

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np, t_eval=t)

Upvotes: 3

Related Questions