Using solve_ivp inside root or fsolve gives ValueError: setting an array element with a sequence

Question

I need to solve a system of equations with a certain parameter p, and then I need to find the value of p that gives me the desired results. My code looks like (in a simplified version)

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import root

def system(t, y, alpha):
    phi, psi, N = y
    dphi_dt = psi
    dpsi_dt = -3 * H(phi, psi, alpha) * psi - dV_dphi(phi, alpha)
    dN_dt = H(phi, psi, alpha)
    return [dphi_dt, dpsi_dt, dN_dt]

# initial conditions
psi0 = -dV_dphi(phi0, alpha)/(3*H(phi0,0, alpha))  # Initial value of psi
N0 = 0.0    # Initial value of N

# time span
t0 = 0
dt = 0.01  # Time step
t = t0

# Initialize lists to store the results
t_values = [t0]
phi_values = [phi0]
psi_values = [psi0]
N_values = [N0]


# Initial conditions
y = [phi0, psi0, N0]

# Solve the system of differential equations for one step
sol = solve_ivp(system, [t, t+dt], y, args=(alpha,))
y = sol.y[:, -1]
t = sol.t[-1]

# Extract solutions
phi, psi, N = y

# calculate extra info
# save data into a dictionary (containing lists) e.g.
data = {'H':[1,2,3,4,5]}
return data

This code works perfectly and gives me the desired results.

Later on, I need to fix the parameter alpha and find what phi0 gives me data['H'][0] == h, for example. I'm using

def H(phi0):
    data = inflate(phi0, alpha)
    p = data['H'][0] - h
    return p

sol = root(H, 10.0)

where data['H'] is one of the lists I calculate before. The function H returns a number, as expected. When calling for root (or fsolve, the result is the same), the code calls for inflate, which then calls for solve_ivp. At this moment, solve_ivp complains about:

def check_arguments(fun, y0, support_complex):
      5     """Helper function for checking arguments common to all solvers."""
----> 6     y0 = np.asarray(y0)

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

The call for the inflate function never even reaches the return point. solve_ivp is not concluded even once. From what I understand, somehow inside the iterations of solve_ivp a sequence is being passed as an argument inside the initial conditions y0. This only happens when solve_ivp is being called inside root, and not if I call solve_ivp separately.

I have no idea why this is happening or how to solve this. Any help is appreciated

I tried changing root for fsolve, and solve_ivp for odeint. The same error was returned, setting array element with a sequence. but the source of the error was different. I've also tried verifying that the function H works properly and returns a number, which indeed it does. From reading the documentation, I'm sure the error comes when solve_ivp is called inside inflate. The function inflate never reaches its returning point, so only one iteration is being made inside root. I'm sure that solve_ivp works properly inside inflate, since it does give the desired results. The problem happens only when it is called inside root.

I tried changing the solve_ivp method, but some methods gave me the same error, while others took so damn long to finish that I gave up.

Márton Horváth · Accepted Answer · 2024-08-11 20:02:48Z

I think, I can see the problem here but I can not provide a complete code since you omitted the definition of the 'H' and 'dV_dphi' functions from the copied script.

So, in your script 'phi0' and 'psi0' are defined as numpy arrays and this is what causes the problem when 'solve_ivp' tries to convert y0 to a numpy array.

To resolve this, it should be ensured that 'N0', 'phi0' and 'psi0' are scalars, not numpy arrays. For this purpose, 'float' can be used.

After defining 'system', I would continue the code as follows, incorporating these changes:

    ###your code including def system(t, y, alpha)...

    def inflate(phi0, alpha):
    phi0 = float(phi0)  # Ensure phi0 is a scalar
    psi0 = float(-dV_dphi(phi0, alpha) / (3 * H(phi0, 0, alpha)))  # Ensure psi0 is a scalar
    N0 = 0.0  # Initial value of N
    
    # initial conditions
    y = [phi0, psi0, N0]

    # time span
    t0 = 0
    dt = 0.01  # Time step
    t = t0


    print(f"Initial conditions: y = {y}")


    sol = solve_ivp(system, [t, t + dt], y, args=(alpha,))
    y = sol.y[:, -1]
    t = sol.t[-1]


    phi, psi, N = y


    data = {'H': [H(phi, psi, alpha)]} 


    print(f"Solved values: phi = {phi}, psi = {psi}, N = {N}")
    print(f"Data: {data}")

    return data

def H_function(phi0, alpha, h):
    data = inflate(phi0, alpha)
    p = data['H'][0] - h
    return p

# Example values 
alpha = 1.0
h = 1.0

# finding the optimal phi0
sol = root(H_function, 10.0, args=(alpha, h))

print(f"Root finding result: {sol}")


# extracting the optimal phi0
optimal_phi0 = sol.x[0]
print(f"Optimal phi0: {optimal_phi0}")

I used arbitrary definitions for the above-mentioned functions and the code worked.

Collectives™ on Stack Overflow

Using solve_ivp inside root or fsolve gives ValueError: setting an array element with a sequence

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related