Why does abs cause bad matlabfunction evaluation?

If I run

x = sym('x','real');
y = sym('y','real');
f = abs(((x - 1)^2 + y^2)^(1/2) - 1)^2 + ((x^2 + y^2)^(1/2) - 1)^2;
d2fdx2 = diff(diff(f,x),x);
[eval(subs(subs(d2fdx2,x,1),y,1)) feval(matlabFunction(d2fdx2),1,1)]

Then I see

ans =

    1.2929       NaN

If I remove the abs then I get the correct evaluation for both ([1.2929 1.2929]).

I'm wondering why this abs is causing trouble in matlabFunction . (in this minimal example it's of course easy to remove)

[Alec Jacobson]

Symbolically we know that if x and y are real then ((x - 1)^2 + y^2)^(1/2) is real, but perhaps the symbolic toolbox is not acknowledging this?

[Dyuman Joshi]

It is acknowledging it -

syms x y real
assumptions

ans = 

eqn = ((x - 1)^2 + y^2)^(1/2) 

eqn = 

isAlways(eqn>=0)

ans = logical
   1

Answer 1

[Catalytic]

The real problem appears to be that f is not fully simplified, and therefore the Symbolic Toolbox tries to naively apply the Chain rule to sub-expressions abs(z)^2. Even though abs(z)^2=z^2 is differentiable for real z, the chain rule isn't applicable because the nested function abs(z) is not continuously differentiable.

This causes the Toolbox to express the derivatives using Dirac deltas which evaluate to NaN at x=y=1 even though that is fake calculus -

x = sym('x','real');
y = sym('y','real');

f = abs(((x - 1)^2 + y^2)^(1/2) - 1)^2 + ((x^2 + y^2)^(1/2) - 1)^2

f = 

d2fdx2 = diff(diff(f,x),x)  %Dirac delta in 6th term.

d2fdx2 = 

On the other hand, when f is properly simplified, the Chain rule works properly, because there are no abs(z)^2 sub-expressions -

fsimp=simplify(  abs(((x - 1)^2 + y^2)^(1/2) - 1)^2 ) + ((x^2 + y^2)^(1/2) - 1)^2

fsimp = 

d2fdx2 = diff(diff(fsimp,x),x)  %no Diracs

d2fdx2 = 

eval(subs(d2fdx2,[x,y],[1,1]))

ans = 1.2929

Answer 2

[Walter Roberson]

You are taking the derivative of abs(), but abs() does not have a continuous derivative.

You have the sub-expression abs(((x - 1)^2 + y^2)^(1/2) - 1)^2 but as x and y both approach 1, then that goes to abs(0) and the derivative of abs at 0 is not continuous.

Depending on the order of substitution, as you partially evaluate the derivative, the Symbolic Toolbox might turn part of the expression into a dirac delta. But when that dirac delta is evaluated at 0 then the result is infinity; that infinity happens to be multiplied by a 0 in that case, and that leads to 0*inf which is NaN.

[eval(subs(subs(d2fdx2,x,1),y,1)) feval(matlabFunction(d2fdx2),1,1)]

There is no documented meaning for eval() of a symbolic expression. You will not find eval() of a symbolic expression described anywhere in doc .

Symbolic expressions are represented to users in a human-readable form rather than as MATLAB expressions or as expressions in the internal symbolic computing language MuPAD. eval() of them can fail.

Answer 3

[Matt J]

It is not the matlabFunction version that is evaluating things badly. It is the first version, with the double application of subs. The comparison that you should be doing is,

x = sym('x','real');
y = sym('y','real');
f = abs(((x - 1)^2 + y^2)^(1/2) - 1)^2 + ((x^2 + y^2)^(1/2) - 1)^2;
d2fdx2 = diff(diff(f,x),x);


[eval(subs(d2fdx2,[x,y],[1,1])) feval(matlabFunction(d2fdx2),1,1)]

ans = 1×2

NaN NaN

which is the correct result.

[Alec Jacobson]

Oh, I didn't know it was an issue to apply `subs` twice. Why is that?

[Matt J]

When subs() is used, a ggeneric simplification is triggered, which is not valid in all cases. Consider the simpler example,

syms x y

expr=(x-1)/(y-1)

expr = 

This expression should be NaN at x=y=1. However, if you first evaluate at x=1, giving the expression 0/(1-y), the symbolic engine simplifies this to its generic value, which is 0,

step1=subs(expr,x,1)

step1 = 

0

and so the second substitution step of course results in 0 as well.

step2=subs(step1,y,1)

step2 = 

0

Conversely,

subs(expr,[x,y],[1,1])

Error using sym/subs
Division by zero.