Determining the roots, maxima, and minima of a discontinuous symbolic function

How can the roots of the following discontinuity function be determined?

syms x 
f =(x^3*heaviside(x))/18 - (3*x^2*heaviside(x))/2 + (27*heaviside(x - 6)*(x - 6))/4 + (27*x*heaviside(x))/4 + heaviside(x - 9)*(x - 9)^2 - (heaviside(x - 9)*(x - 9)^3)/18
fplot(f,[-2,12])

I have tried using feval and the roots functions but neither of these methods worked. I would appreciate any help.

Answer 1

[Shubham Rawat]

Hi Aleem,

Here is a function in syms called solve, which can evaluate roots of a function but for continous part only,

syms x
f =(x^3*heaviside(x))/18 - (3*x^2*heaviside(x))/2 + (27*heaviside(x - 6)*(x - 6))/4 + (27*x*heaviside(x))/4 + heaviside(x - 9)*(x - 9)^2 - (heaviside(x - 9)*(x - 9)^3)/18

solve(f == 0, x, 'MaxDegree', 4); % for solving the equation 
roots = vpa(ans,6);

fplot(f,[-2, 12])
hold on
plot(roots, subs(f,roots), '*')
hold off

For Maxima and Minima,

syms x
f =(x^3*heaviside(x))/18 - (3*x^2*heaviside(x))/2 + (27*heaviside(x - 6)*(x - 6))/4 + (27*x*heaviside(x))/4 + heaviside(x - 9)*(x - 9)^2 - (heaviside(x - 9)*(x - 9)^3)/18

g = diff(f,x); % for differentiation
solve(g == 0,x,'MaxDegree',4); %for solving g

extrema = vpa(ans,6)
fplot(f,[-2, 12])
hold on
plot(extrema, subs(f,extrema), '*')
hold off

This code works fine for the differentiable part, for non-differentiable part you may iterate the points one-by-one how function is behaving.

[Walter Roberson]

Caution: solve(f) returns 9, -1, and 27/2 - (9*3^(1/2))/2 . The 9 and 27/2 - (9*3^(1/2))/2 are exact, but the -1 is "representative".

The function is 0 for x <= 0, so there are an infinite number of roots. When you use solve() in the form shown above, when there is an interval, MATLAB displays a "representative" value from the interval. In the past I have posted a relatively detailed list of how it chooses the representative value, but the details do not really matter at this point.

In order to get a symbolic representation of the interval, you would need

sol = solve(f, 'returnconditions', true);
sol.x
sol.conditions

and you will see that one of the solutions is x and that the corresponding condition is x <= 0

[Aleem Andrew]

Thanks for your help