Problem with _plus() and symbolic differentiation

I define a symbolic expression (dx_3dt) and make an assumption that reduces that expression to one term. By this the function _plus(...) appears in my expression (Output dx_3dt) which is apparently not differentiable by Matlabs definition leading to NaN value. Differentiating the same expression "manually" allows me to obtain the correct solution.

Any advice how to handle this?

MWE:

clear all;
reset(symengine)

theta = sym('theta_%d', [5, 1], {'positive', 'real'}); % physical parameters of system dynamics
x = sym('x_%d', [6, 1], {'positive', 'real'}); % system states
grav = 9810;  % gravitation accelaration [mm/s^2]

eps = 1e-10;
outflow = piecewise(x(3) <= eps, 0, -theta(5)*sqrt(2*grav*x(3)));
dx_3dt = theta(4)*sign(x(2)-x(3))*sqrt(2*grav*abs(x(2)-x(3)))+theta(2)*x(5)+outflow;

simplify(diff(dx_3dt, x(2)))
assumeAlso(x(2) == x(3))
simplify(dx_3dt)

Output dx_3dt in command window:

piecewise(x_3 <= 1/10000000000, _plus(theta_2*x_5), theta_2*x_5 - 6*545^(1/2)*theta_5*x_3^(1/2) + 6*545^(1/2)*theta_4*abs(x_2 - x_3)^(1/2)*sign(x_2 - x_3))

Expected result:

simplify(diff(dx_3dt, x(2)))
diff(theta(2)*x(5), x(2))

First ideas

The problem seems to be the piecewise function since omitting piecewise gives the expected result:

clear all;
reset(symengine)

theta = sym('theta_%d', [5, 1]); % physical parameters of system dynamics
x = sym('x_%d', [6, 1]); % system states
grav = 9810;  % gravitation accelaration [mm/s^2]

eps = 1e-10;
% outflow = piecewise(x(3) <= eps, 0, -theta(5)*sqrt(2*grav*x(3)));
dx_3dt = theta(4)*sign(x(2)-x(3))*sqrt(2*grav*abs(x(2)-x(3)))+theta(2)*x(5)-theta(5)*sqrt(2*grav*x(3));

simplify(diff(dx_3dt, x(2)))
assumeAlso(x(2) == x(3))
simplify(dx_3dt)

simplify(diff(dx_3dt, x(2)))
diff(theta(2)*x(5), x(2))

Problem with diff and assumptions

Furhtermore if I apply the following workaround with replacing the critical term with a symbolic function of x_3 and assume that x_3 == x_2 and differentiate with respect to x_2 it is zero! So the diff function appears to also not consider assumptions here:

clear all;
reset(symengine)

theta = sym('theta_%d', [5, 1], {'positive', 'real'}); % physical parameters of system dynamics
x = sym('x_%d', [6, 1], {'positive', 'real'}); % system states
grav = 9810;  % gravitation accelaration [mm/s^2]


eps = 1e-10;
% outflow = symfun(x_3, x_3)
syms x_3 outflow(x_3)
% outflow = piecewise(x(3) <= eps, 0, -theta(5)*sqrt(2*grav*x(3)));
dx_2dt = theta(3)*x(6)*sign(x(1) - x(2))*sqrt(2*grav*abs(x(1)-x(2)))-theta(4)*sign(x(2)-x(3))*sqrt(2*grav*abs(x(2)-x(3)))
dx_3dt = theta(4)*sign(x(2)-x(3))*sqrt(2*grav*abs(x(2)-x(3)))+theta(2)*x(5)+outflow(x_3)
f = [dx_2dt; dx_3dt];

simplify(diff(dx_3dt, x(2)))
simplify(diff(dx_3dt, x(3)))
assumeAlso(x(2) == x(3))
simplify(dx_3dt)

assumptions(x(2))
assumptions(x(3))
simplify(dx_3dt)
simplify(diff(dx_3dt, x(2)))

simplify(diff(dx_3dt, x(3)))

jac = simplify(jacobian(f, x))

Answer 1

[Stefanie Schwarz]

Thank you for bringing this to our attention. I have alerted our developers to this issue, as of course, there should appear no "_plus" in the output.

The workaround is, as you noticed, to use "simplify" or any other operation that triggers a re-evaluation, such as dx_3dt * 1.

Please note that an assumption x(2) == x(3) does not make x(2) "depend" on x(3), and it does not have the effect of an assignment . I.e., even if you assume x == 42, the derivative of x^2 is 2*x while the derivative of 42^2 is 0 and the derivative diff(42^2, 42) makes no sense. This is similar for integrals and other cases of bound variables (summation variables etc.).