CODEWITHSUNDEEP
pythonscipyodeint
mg547
mg547

Reputation: 55

Using vector inputs to odeint in python scipy to solve a system of two differential equations

I am trying to solve a system of two coupled differential equations using python odeint().

In this system, a function f depends on two variables f(y,t) and another function g depends on one variable g(t). For example, something like that (illustrative only):

I have tried using the following code:

import numpy as np
from scipy import integrate

dy = 0.05
y = np.arange(0, 1 + dy, dy)
dt = 1
tmax = 100

t = np.arange(0,tmax,dt)
f = np.ones([tmax,len(y)])
g = np.ones(tmax)

def deriv(y,t):
    fi = y[0]
    gi = y[1]

    fprime = (1 - y) + fi + gi
    gprime = gi

    return [fprime, gprime]

# Initial conditions
f_ini = np.ones(len(y))*15          
g_ini = np.array([0.3])
sol_ini = np.concatenate((f_ini, g_ini), axis=0)

# solve the DEs
soln = integrate.odeint(deriv, sol_ini, t)

I get the following error at the last line of my code:

ValueError: setting an array element with a sequence.

I am guessing that I am not setting my initial conditions correctly. Any advice?

Upvotes: 2

Views: 2375

Answers (1)

Dietrich
Dietrich

Reputation: 5531

You might not have chosen the best set of equations: Your equations are independent of each other and the first one is a partial differential equation. With Sympy you're able to find closed form solutions:

from IPython.display import display
import sympy as sy
from sympy.solvers.ode import dsolve
from sympy.solvers.pde import pdsolve


sy.init_printing()  # LaTeX like pretty printing for IPython


t, y = sy.symbols("t, y", real=True)
f, g = sy.symbols("f, g", function=True)

eq1 = sy.Eq(g(t).diff(t), g(t))
g_sol = dsolve(eq1)
print("For the ode")
display(eq1)
print("the solution is")
display(g_sol)

eq2 = sy.Eq(f(y, t).diff(t), (1-y) + f(y, t) + g_sol.rhs)
f_sol = pdsolve(eq2)
print("For the pde")
display(eq2)
print("the solution is")
display(f_sol)

gives in an IPython interpreter

solution

You see that g(t) has an undetermined constant C_1 and f(t) and undetermined function F(y), which need to be defined by initial conditions. If you know what your system behaves like at, e.g. time t=0, C_1 and F(y) can easily be determined.

Upvotes: 1

Related Questions