Solving a 2-Variable function to only 1 of the variable doesnt work on arrays

Question

I have found from StackOverflow how to solve a fucntion with two variables giving the one as constant/known.

This is the part of the code:

def R(gg,a):
    return a-r0*g0**(1/2)*D(gg)/gg**(1/2)
def G(r):
    partial_func = functools.partial(R, a=r)
    return fsolve(partial_func,10,xtol=10**-1)

and it works, since for the first 2 prints, I get the same value

f=([10,15])
print(G(10))
print(G(f[0]))
print(G(f))

but when giving the full array it has the following error:

The array returned by a function changed size between calls

Are you aware that `f` is a tuple of a list?? I think you mean `f=[10, 15]`..right? — Anwarvic, Jan 27 '19 at 13:10
@Anwarvic It's actually a list since a one element tuple is `([10,15],)` but the brackets are not necessary. — Jacques Gaudin, Jan 27 '19 at 13:12
`f` is not a tuple of a list. it would be if you would write `f=([10,15], )` — Jonathan R, Jan 27 '19 at 13:12
@billy You should attach the returned values and definitions for all variables to make it easier to help — Jonathan R, Jan 27 '19 at 13:19
Per [the docs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fsolve.html#scipy.optimize.fsolve), `fsolve` expects `partial_func` to take a vector as input and return a vector of the same length. It looks like `partial_func` is taking a scalar but returning a vector. — unutbu, Jan 27 '19 at 13:19

Jacques Gaudin · Accepted Answer · 2019-01-27T14:01:26.377

2

It looks like you are trying to find the roots of R for different values included in f.

The problem is that partial_func has an single value as starting estimate and wants to return an array of the same length as a (in your case 2 values).

In other words, there is not a single value root to your problem. For example the root for f[0]=10 is probably different from the root for f[1]=15. The solution should be an array of two values in this case.

To fix this, you need to give an array for the x0 (starting estimate) parameter of fsolve.

def G(r):
    partial_func = functools.partial(R, a=r)
    return fsolve(partial_func,[10,10],xtol=10**-1)

So that for each values in a there is a initializer for gg and the solution is a vector of the same length as f.

edited Jan 27 '19 at 14:01

answered Jan 27 '19 at 13:39

Jacques Gaudin

15,779
10
54
75

Wait, actually I have a problem. When I give G(f) as a result I get ([f[0],f[1]]) and not G(f[0]),G(f[1]). What's to blame here? :( – Billy Matlock Jan 27 '19 at 15:41
I don't understand. See [this snippet](https://codepad.co/snippet/1IkWsM9g) which returns the roots of `R` for the values present in `f`. Let me know. – Jacques Gaudin Jan 27 '19 at 16:15
I get `f=[10 15 17] G(f[0])= [30.73218835] G(f[1])=[39.49923357] G(f[2])=[41.79696086]` and `G(f)=[30.73259241 39.66454009 42.16416723]~=[G(f[0], G(f[1], G(f[2])]` – Jacques Gaudin Jan 27 '19 at 16:34
Any reason why they are not the same values? – Billy Matlock Jan 27 '19 at 16:40
There is no guarantee that you will have the exact same values. This is a numerical approach where "The calculation will terminate if the relative error between two consecutive iterates is at most xtol." If you refine `xtol` the algorithm will run longer and give you more precise and similar answers. – Jacques Gaudin Jan 27 '19 at 17:04
Just a last code to see that there is a problem. If you run [this](https://pastebin.com/VbcmPHM8) you will see that there is a huge difference. I mean it can't be the tolerance. :( :( – Billy Matlock Jan 27 '19 at 17:19
There is indeed a huge difference because the `e00` function is very flat and close to zero. That's the limit of a numerical approach. – Jacques Gaudin Jan 27 '19 at 17:43
BTW, carefully choosing the starting estimate is very important. You will see that there is an error message in the result of your last example if you pass `full_output=True` – Jacques Gaudin Jan 27 '19 at 18:09

score 1 · Answer 2 · answered Jan 27 '19 at 13:21

So reading this w/o knowing all parameters used in the function i would say that in case of print(G(f)) you provide a scalar and return an array, which does not work.

Try calling your function with a=f and look at the returned value.

The docs state:

fsolve: func: A function that takes at least one (possibly vector) argument, and returns a value of the same length

Solving a 2-Variable function to only 1 of the variable doesnt work on arrays

2 Answers2