MATLAB: Find the root of a double integral

root of the 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)

Best Answer

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.
Related Question