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

-1

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.

edited Jan 3 at 20:31

asked Jan 2 at 7:06

kemal kaya

1

When you ask about an error, you need to including as much information as possible - the full error message (not just the title) plus the traceback (the location information). We can't guess where the problem occurs (without running the code outselves).

– hpaulj
Jan 2 at 7:16

@hpaulj i've added the error message.

– kemal kaya
Jan 2 at 7:37

It's hard to trace all that's going on, but it appears that when evaluating jac_np a calculation produces 3 values, but it expects 4 (this unpacking business). I suppose the 3 comes from U_0. I would try to evaluate jac_np along, outside of the ode calc, e.g. jac_np(0, U0).

– hpaulj
Jan 2 at 17:22

@hpaulj i've simplified my code. As you suggested I evaluated jac_np along outside the ode calc and seems no problem. but sending array arguments to jac_np causes error.

– kemal kaya
Jan 3 at 20:36

add a comment |

-1

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.

edited Jan 3 at 20:31

asked Jan 2 at 7:06

kemal kaya

1

When you ask about an error, you need to including as much information as possible - the full error message (not just the title) plus the traceback (the location information). We can't guess where the problem occurs (without running the code outselves).

– hpaulj
Jan 2 at 7:16

@hpaulj i've added the error message.

– kemal kaya
Jan 2 at 7:37

It's hard to trace all that's going on, but it appears that when evaluating jac_np a calculation produces 3 values, but it expects 4 (this unpacking business). I suppose the 3 comes from U_0. I would try to evaluate jac_np along, outside of the ode calc, e.g. jac_np(0, U0).

– hpaulj
Jan 2 at 17:22

@hpaulj i've simplified my code. As you suggested I evaluated jac_np along outside the ode calc and seems no problem. but sending array arguments to jac_np causes error.

– kemal kaya
Jan 3 at 20:36

add a comment |

-1

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.

edited Jan 3 at 20:31

asked Jan 2 at 7:06

kemal kaya

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.

scipy sympy ode

edited Jan 3 at 20:31

asked Jan 2 at 7:06

kemal kaya

edited Jan 3 at 20:31

asked Jan 2 at 7:06

kemal kaya

edited Jan 3 at 20:31

asked Jan 2 at 7:06

kemal kaya

asked Jan 2 at 7:06

kemal kaya

asked Jan 2 at 7:06

kemal kaya

1

When you ask about an error, you need to including as much information as possible - the full error message (not just the title) plus the traceback (the location information). We can't guess where the problem occurs (without running the code outselves).

– hpaulj
Jan 2 at 7:16

@hpaulj i've added the error message.

– kemal kaya
Jan 2 at 7:37

It's hard to trace all that's going on, but it appears that when evaluating jac_np a calculation produces 3 values, but it expects 4 (this unpacking business). I suppose the 3 comes from U_0. I would try to evaluate jac_np along, outside of the ode calc, e.g. jac_np(0, U0).

– hpaulj
Jan 2 at 17:22

@hpaulj i've simplified my code. As you suggested I evaluated jac_np along outside the ode calc and seems no problem. but sending array arguments to jac_np causes error.

– kemal kaya
Jan 3 at 20:36

add a comment |

1

When you ask about an error, you need to including as much information as possible - the full error message (not just the title) plus the traceback (the location information). We can't guess where the problem occurs (without running the code outselves).

– hpaulj
Jan 2 at 7:16

@hpaulj i've added the error message.

– kemal kaya
Jan 2 at 7:37

It's hard to trace all that's going on, but it appears that when evaluating jac_np a calculation produces 3 values, but it expects 4 (this unpacking business). I suppose the 3 comes from U_0. I would try to evaluate jac_np along, outside of the ode calc, e.g. jac_np(0, U0).

– hpaulj
Jan 2 at 17:22

@hpaulj i've simplified my code. As you suggested I evaluated jac_np along outside the ode calc and seems no problem. but sending array arguments to jac_np causes error.

– kemal kaya
Jan 3 at 20:36

When you ask about an error, you need to including as much information as possible - the full error message (not just the title) plus the traceback (the location information). We can't guess where the problem occurs (without running the code outselves).

– hpaulj
Jan 2 at 7:16

@hpaulj i've added the error message.

– kemal kaya
Jan 2 at 7:37

It's hard to trace all that's going on, but it appears that when evaluating jac_np a calculation produces 3 values, but it expects 4 (this unpacking business). I suppose the 3 comes from U_0. I would try to evaluate jac_np along, outside of the ode calc, e.g. jac_np(0, U0).

– hpaulj
Jan 2 at 17:22

@hpaulj i've simplified my code. As you suggested I evaluated jac_np along outside the ode calc and seems no problem. but sending array arguments to jac_np causes error.

– kemal kaya
Jan 3 at 20:36

add a comment |

1 Answer
1

active

oldest

votes

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)

answered Jan 3 at 22:02

LutzL

14.2k21426

add a comment |

Your Answer

StackExchange.ifUsing("editor", function () {
StackExchange.using("externalEditor", function () {
StackExchange.using("snippets", function () {
StackExchange.snippets.init();
});
});
}, "code-snippets");

StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "1"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);

StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});

function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});

}
});

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fstackoverflow.com%2fquestions%2f54002500%2fhow-to-use-a-sympy-generated-jacobi-matrix-in-the-solution-of-an-ode-system%23new-answer', 'question_page');
}
);

Post as a guest

Name

Required, but never shown

1 Answer
1

active

oldest

votes

1 Answer
1

active

oldest

votes

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)

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)

answered Jan 3 at 22:02

LutzL

14.2k21426

add a comment |

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)

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)

answered Jan 3 at 22:02

LutzL

14.2k21426

add a comment |

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)

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)

answered Jan 3 at 22:02

LutzL

14.2k21426

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)

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)

answered Jan 3 at 22:02

LutzL

14.2k21426

answered Jan 3 at 22:02

LutzL

14.2k21426

answered Jan 3 at 22:02

LutzL

14.2k21426

answered Jan 3 at 22:02

LutzL

14.2k21426

add a comment |

draft saved

draft discarded

Thanks for contributing an answer to Stack Overflow!

Please be sure to answer the question. Provide details and share your research!

But avoid …

Asking for help, clarification, or responding to other answers.

Making statements based on opinion; back them up with references or personal experience.

To learn more, see our tips on writing great answers.

draft saved

draft discarded

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Sign up or log in

StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});

Post as a guest

Name

Required, but never shown

Name

Required, but never shown

Name

Required, but never shown

This page is only for reference, If you need detailed information, please check here

Search This Blog

Ufyukyu