Find the root of a double integral

why the following code cannot find the root of the function g(z). Note f(2)=6. So fzero(g,1) should return 2.

syms x y z
f(z) = vpaintegral(vpaintegral(x*y, x, [1 3]), y, [-1 z]);
f(2)
g(z) = f(z) - f(2);
fzero(g,1)

Answer 1

[Walter Roberson]

fzero is not defined for symbolic functions -- but fzero also does not test for that case.

The code sees that symbolic functions have an eval() method and thinks that everything must be okay, and starts processing. It evaluates at two locations, and then tries

if (fa > 0) ~= (fb > 0)

At that point in the code, fa and fb are symbolic numbers that are both negative. But in MATLAB, (symbolic RELATIONSHIP symbolic) does not evaluate to boolean, even when the symbolics are both symbolic numbers. So internally this ends up testing

if (0.0 < -6.1115370849898476459415014861013) ~= (0.0 < -6.0)

and since the two are not the same expression on both sides, the ~= is deemed to hold, and it thinks it has found a sign change. A couple more rounds in which it keeps thinking it has found a sign change, until the interval is narrow enough, and out pops 0.971715728752539

The solution is to not use fzero for this purpose. Use vpasolve() instead, as vpasolve() is able to accept the range hint that you are passing. solve() is also able to find a solution, but there are two solutions -2 and +2 and since you cannot pass the range hint as such, solve() happens to find the negative solution.

The alternative is

syms z positive
solve(g)

This gives the needed range hint.

[Yufei Cao]

Dear Walter,

The solution works very well! Thanks for the detail explanation!