Minimizing function in terms of integral - scipy

Question

I want to find the parameter of a function that will result in a specific integral value in a defined interval. To make things simpler, the example below considers the function to be a straight line, and the parameter I want to find is the slope m. I use scipy.integrate.quad to integrate and was trying to use scipy.integrate.minimize to find the slope:

from scipy.integrate import quad
from scipy.optimize import minimize

IntegralValue = 10
initial_guess = 0.1

def function(m):
    intgrl, abserr = quad(lambda x: m*x, 0, 10) # interval from 0 to 10
    return intgrl - IntegralValue
res = minimize(function, initial_guess, method='nelder-mead',options={'xtol': 1e-8, 'disp': True})

print(res.x)

A simple 0.2 should be the result (function(0.2) results in 0), but for some reason I get a warning that the "maximum number of function evaluations has been exceeded" and the output is something that doesn't make sense (-6.338253e+27). Default minimize method also did not work. Maybe I am missing something or misunderstood, but I couldn't figure it out. I appreciate the help. I am using Python 3.7.10 on Windows 10.

Dumb question: are you sure that zero is the _minimum_ of the `function`? Can `function` return negative values? I think `intgrl - IntegralValue` in general can easily be negative. If it can, `minimize` will try to find the lowest _negative_ value of this function, which may be negative infinity. — ForceBru, Mar 07 '22 at 11:27
If that's a dumb question, I am even dumber. It makes total sense what you said and I think that is the problem... that means I have to put a constraint for non-negative value somehow, somewhere, right? Maybe just an ```abs()```... — andrerud, Mar 07 '22 at 11:30
Hint: minimizing `(intgrl - IntegralValue)**2` is equivalent to solving the equation `integrl = integralValue`. — joni, Mar 07 '22 at 11:31
I think the simplest thing is to return a different kind of distance between the two values: `abs(intgrl - IntegralValue)` or `(intgrl - IntegralValue)**2` maybe. — ForceBru, Mar 07 '22 at 11:31
Thank you both, yea, I knew it was going to be something obvious — andrerud, Mar 07 '22 at 11:34

andrerud · Accepted Answer · 2022-03-07T12:37:10.577

1

The minimum is of course negative infinity. To get a zero, I just had to constrain the problem to non negative values, as commented by @ForceBru and @joni. Thanks!

EDIT

here is the final functional full code:

from scipy.integrate import quad
from scipy.optimize import minimize

IntegralValue = 10
initial_guess = 0.1

def function(m):
    intgrl, abserr = quad(lambda x: m*x, 0, 10) # interval from 0 to 10
    return abs(intgrl - IntegralValue)
res = minimize(function, initial_guess, method='nelder-mead',options={'xtol': 1e-8, 'disp': True})

print(res.x)

edited Mar 07 '22 at 12:37

answered Mar 07 '22 at 11:38

andrerud

146
1
11

2

It may be helpful for future generations to include a snippet of code that shows how to fix the issue - maybe a definition of `function` which now returns `abs` or `(...)**2` – ForceBru Mar 07 '22 at 11:41

score 1 · Answer 2 · answered Mar 07 '22 at 14:26

The problem that you describe is one of root finding; that is, you want to find m such that function(m) = 0. The appropriate tool for such a problem is one of the root finding functions in SciPy.

For example, you could use root_scalar as in the following code, which uses your original function function:

from scipy.integrate import quad
from scipy.optimize import root_scalar


def function(m):
    intgrl, abserr = quad(lambda x: m*x, 0, 10) # interval from 0 to 10
    return intgrl - IntegralValue


IntegralValue = 10
guess_bracket = [0.0, 1.0]

res = root_scalar(function, x0=guess_bracket[0], x1=guess_bracket[1])

if res.converged:
    print(f"m = {res.root}")
else:
    print(f'root_scalar failed: res.flag={res.flag}')

Output:

m = 0.2

Minimizing function in terms of integral - scipy

2 Answers2