Why does abs cause bad matlabfunction evaluation?

3 次查看(过去 30 天)
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)
  2 个评论
Alec Jacobson
Alec Jacobson 2023-12-9
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
Dyuman Joshi 2023-12-9
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

请先登录,再进行评论。

回答(3 个)

Catalytic
Catalytic 2023-12-9
编辑:Catalytic 2023-12-9
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

Walter Roberson
Walter Roberson 2023-12-9
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.

Matt J
Matt J 2023-12-9
编辑:Matt J 2023-12-9
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.
  2 个评论
Alec Jacobson
Alec Jacobson 2023-12-9
Oh, I didn't know it was an issue to apply `subs` twice. Why is that?
Matt J
Matt J 2023-12-9
编辑:Matt J 2023-12-9
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.

请先登录,再进行评论。

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by