How to solve a system of ODEs with Scipy when the variables in the equations are autogenerated-CodePudding

I'm generating a system of ODEs in the form of a list of equations where the variables are sympy.Symbol(), for example [3*sympy.Symbol('x')**3 sympy.Symbol('y') , sympy.Symbol('x')-sympy.Symbol('y')**4].

So using this example, I can solve this system through the code

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

def kin(t, z):
    x, y = z
    eqs = [3*x**3 y  , x-y**4]
    return eqs

y0 = [0, 10]
t = np.linspace(0, 10000, 10000)
sol = solve_ivp(kin, [0, 10000], y0, dense_output=True)
print(sol.sol(t).T)
plt.plot(t, sol.sol(t).T)
plt.show()

And it will give this result https://imgur.com/3KvOkrU

But since the variables of x and y are in reality generated by a different code, I figured I could do

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

dic = {'x': Symbol('x'), 'y': Symbol('y')}
vars = ('x', 'y')
#eqs =[3x^3 y , x-y^4]
eqs = [3*dic[vars[0]]**3 dic[vars[1]], dic[vars[0]]-dic[vars[1]]**4]

def kin(t, z):
    dics = {}
    for v in range(len(vars)):
        dics[vars[v]] = z[v]

    for e in range(len(eqs)):
        eqs[e] = eqs[e].subs(dics)
    return eqs

y0 = [0, 10]
t = np.linspace(0, 10000, 10000)
sol = solve_ivp(kin, [0, 10000], y0, dense_output=True)
print(sol.sol(t).T)
plt.plot(t, sol.sol(t).T)
plt.show()

Where it merely replaces the values of x and y from the first part by z[0] and z[1] but the output of this code is completely different https://imgur.com/ovi3mk8

I dont see where the issue is and would appreciate if anyone code point out where the problem is or help me find a way to evaluate these auto-generated systems of ODEs .

CodePudding user response：

The reason you are getting that results is because you are overwriting eqs inside kin at each iteration. Also, as mentioned by Lutz, you should use lambdify which is going to convert symbolic expressions to numerical functions so that they can be evaluated by Numpy (much faster).

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

dic = {'x': Symbol('x'), 'y': Symbol('y')}
vars = ('x', 'y')
eqs = [3*dic[vars[0]]**3 dic[vars[1]], dic[vars[0]]-dic[vars[1]]**4]
# convert symbolic expressions to numerical functions
eqs = [lambdify(list(dic.values()), e) for e in eqs]

def kin(t, z):
    # evaluate each numerical function with the current values
    return [e(*z) for e in eqs]

y0 = [0, 10]
t = np.linspace(0, 10000, 10000)
sol = solve_ivp(kin, [0, 10000], y0, dense_output=True)
plt.plot(t, sol.sol(t).T)
plt.show()

This will give you the first picture you posted.