ValueError: `fun` must return at most 1-d array_like. f0.shape: (2, 100)-CodePudding

I'm trying to find the intersection of lines after solving a system of equations, and it shows 'ValueError: fun must return at most 1-d array_like. f0.shape: (2, 100)'. Not really sure where I got wrong.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import least_squares
import scipy.interpolate, scipy.optimize


k1 = 10 ** (-9.25)
k2 = 10 ** (-18.26)
k3 = 10 ** (-35.91)
logk = 11.96
pH_x = np.linspace(0, 14, 100)
ni_total = 0.1
n_total = 1.2
T_ = 298
h = 10 ** (-pH_x)
def equations(p):
    nh4, nio2 = p
    f = 10 ** (logk - 2 * pH_x)
    ni2pfree = f/(1   f / ni_total)
    nh3 = k1 * nh4 / h
    nin4p2 = k2 * (nh4 ** 4) / (h ** 2)
    nin6p2 = k3 * (nh4 ** 6) / (h ** 4)
    return (n_total-nh3-4*nin4p2-6*nin6p2-nh4,
            ni_total - ni2pfree - nin4p2- nin6p2- nio2)
res = least_squares(equations, (0.1, 0.1), bounds=((0, 0), (n_total,ni_total )),method='dogbox',xtol=1e-12)
nh4 =res.x[0]
nio2 =res.x[1]
f = 10 ** (logk - 2 * pH_x)
ni2pfree = f / (1   f / ni_total)
nh3 = k1 * nh4 / h
nin4p2 = k2 * (nh4 ** 4) / (h ** 2)
nin6p2 = k3 * (nh4 ** 6) / (h ** 4)

y1 = -0.2405   0.0296 * np.log10(ni2pfree)
y2 = 0.1102-0.0592*pH_x
interp1 = scipy.interpolate.InterpolatedUnivariateSpline(pH_x, y1)
interp2 = scipy.interpolate.InterpolatedUnivariateSpline(pH_x, y2)

def difference(pH_x):
    return np.abs(interp1(pH_x) - interp2(pH_x))

x_at_crossing = scipy.optimize.fsolve(difference, x0=0.0)
plt.plot(x_at_crossing, interp1(x_at_crossing), 'cd', ms=7, label='fsolve method')
plt.show()

CodePudding user response：

The function least_squares() expect the first argument fun with the signature fun(x, *args, **kwargs) where the argument x is an ndarray of shape (n,), while in your case it has shape (2,100).

This error arises from the way you return the function equation():

return (n_total-nh3-4*nin4p2-6*nin6p2-nh4, ni_total - ni2pfree - nin4p2- nin6p2- nio2) returns a tuple of shape (2,) and not an array of shape (n,).

Changing your return line to:

return np.array([n_total-nh3-4*nin4p2-6*nin6p2-nh4, ni_total - ni2pfree - nin4p2- nin6p2- nio2]).flatten()

would allow the code to run entirely, but I believe the answer you'd get is not correct, since this 'solution' only flattens the array that had shape (2,100) to shape (n,).