How to simplify function handles?

Y =

function_handle with value:

@(s) (F(s)-Q(s))/P(s)

As seen, I have three functions F, Q and P and Y is function made by algebraic operations on these three. Is there a way, to get a simplified expression for Y in terms of the function variable (s) ? When I type in Y in command window, I just get the result as shown above.

As a simplified example:

F = @(s) s;
Q = @(s) s^2;
P = @(s) s^3-3*s^2-1;
Y = @(s) (F(s) - Q(s))/P(s);

When i type in Y in command window, what is see is:

Y =

function_handle with value:

@(s) (F(s)-Q(s))/P(s)

What I wish to see is:

Y = (s-s^2)/(s^3-3*s^2-1)

[the cyclist]

It would helpful if you included a small example of what you see, and what you wish to see.

[anuj maheshwari]

Added

Answer 1

[Walter Roberson]

You used a tag of symbolic which implies that you might have the Symbolic Toolbox. If so then

syms s
Ys = Y(s)

You might want to simplify() the result.

[Walter Roberson]

F = @(s) s;
Q = @(s) s.^2;
P = @(s) s.^3-3*s.^2-1;
Y = @(s) (F(s) - Q(s))./P(s);

syms s real
Ys = Y(s)

Ys = 

dY = diff(Ys, s)

dY = 

critical = solve(dY, s)

critical = 

d2Y = diff(dY, s)

d2Y = 

crit2 = subs(d2Y, s, critical)

crit2 = 

definitely_max = isAlways(crit2 < 0, 'unknown', 'false');
definitely_min = isAlways(crit2 > 0, 'unknown', 'true');
not_sure = ~(definitely_max | definitely_min);

max_locations = critical(definitely_max)

max_locations = 

max_values = Y(max_locations)

max_values = 

vpa(max_values)

ans = 

0.42303605210393856972993261326406

min_locations = critical(definitely_min)

min_locations = 

min_values = Y(min_locations)

min_values = 

vpa(min_values)

ans = 

other_locations = critical(not_sure)

 
other_locations =
 
Empty sym: 0-by-1
 

other_values = Y(other_locations)

 
other_values =
 
Empty sym: 0-by-1
 

[anuj maheshwari]

Thanks a lot !!!

[Walter Roberson]

F = @(s) s;
Q = @(s) s.^2;
P = @(s) s.^3-3*s.^2-1;
Y = @(s) (F(s) - Q(s))./P(s);

syms s real
Ys = Y(s)

Ys = 

fplot(Ys, [-10 10])

[N,D] = numden(Ys)

N = 

D = 

location_of_zeros = solve(N,s)

location_of_zeros = 

location_of_poles = solve(D, s, 'maxdegree', 4)

location_of_poles = 

vpa(location_of_poles)

ans = 

3.1038034027355365331649473328289

I guess my calculas analysis needed to be more careful as it missed the +/- infinity of the pole.

[anuj maheshwari]

syms s

Y = Y(s);

y = matlabFunction(ilaplace(Y));

fplot(y, [1e-4,10e-4]);

It throws the error:

Warning: Error updating FunctionLine.

The following error was reported evaluating the function in FunctionLine update: Unrecognized function or variable 's'.

and the plot is empty. Any Suggestions?

[Walter Roberson]

F = @(s) s;
Q = @(s) s.^2;
P = @(s) s.^3-3*s.^2-1;
Y = @(s) (F(s) - Q(s))./P(s);

syms s real
Ys = Y(s)

Ys = 

iY = ilaplace(Ys)

iY = 

matlabFunction(iY)

Error using symengine
Code generation failed due to unexpected object of type 'RootOf'.

Error in sym/matlabFunction>mup2mat (line 432)
res1 = mupadmex('symobj::generateMATLAB',r.s,ano,spa);

Error in sym/matlabFunction>mup2matcell (line 401)
r = mup2mat(c{1},true,sparseMat);

Error in sym/matlabFunction (line 172)
body = mup2matcell(funs, opts.Sparse);

That is, ilaplace() generated a root() [really a RootOf] and matlabFunction() is not able to deal with those.

In theory as this was a polynomial of degree 3, ilaplace could have generated the exact formula instead of using a root() or RootOf() placeholder... but it didn't, and there is no method provided by Mathworks to "fix-up" the expression afterwards.

A few weeks ago, I posted an attempt to do conversion of these kinds of structures. I will try to dig up the code and see if I can get it to work for your case.

[anuj maheshwari]

Oh my actual "simple case" is :

((5946427424506845*s)/2417851639229258349412352 + (100*exp(-s/400000) - 100)/(s*(exp(-s/100000) - 1)) + 707092764722665/144115188075855872)/(s^2/10000000000 + s/5000 + 1)

and things will get more complex from there. If matlab isnt the tool for calculating inverse laplace of these functions, can you recommend something?

[Walter Roberson]

In the below, eiY is the output to pay attention to: it is the symbolic form of the ilaplace expaned out from sym of roots into closed formula.

The plotting after that is investigative, to try to figure out why the fplot() is showing up with discontinuities.

According to the X / YY values, the matlabFunction() version of eiY is producing values with imaginary component roughly 1e-17 there. Zooming in trying to look for visible infinities or whatever, did not help. But taking some of the locations that produce non-zero imaginary part and evaluating them symbolically to 50 digits shows a mix of pure-rule output and very small imaginary parts; the amplitude of the imaginary parts decreases as the number of digits increases.

So what is happening here is that the matlabFunction() version is suffering from round-off error (and even the symbolic version does too). The root() that is the solution of the ilaplace has a conjugate pair, so in theory the positive and negative imaginary components should cancel, but in practice because of floating point round off, they do not.

In this particular case, at least for positive t, you should be able to just take real() of the values.

format long g
digits(50)
F = @(s) s;
Q = @(s) s.^2;
P = @(s) s.^3-3*s.^2-1;
Y = @(s) (F(s) - Q(s))./P(s);
syms s real
Ys = Y(s)

Ys = 

iY = ilaplace(Ys)

iY = 

Nth = @(expr,N) expr(N);
resolve_a_root = @(X6) Nth(solve(children(X6,1),children(X6,2),'maxdegree',4),children(X6,3));
resolve_a_summand = @(X5) mapSymType(X5, 'root', resolve_a_root);
resolve_nth_summand = @(X2,N) resolve_a_summand(subs(children(X2,1),children(X2,2),N));
map_symsum = @(X2) sum(arrayfun(resolve_nth_summand, repmat(X2,1,double(children(X2,4)-children(X2,3)+1)), children(X2,3):children(X2,4)));
expand_symsum = @(X1) mapSymType(X1, 'symsum', map_symsum);

eiY = simplify(expand_symsum(iY))

eiY = 

y = matlabFunction(eiY);
fplot(y, [1e-4,10e-4]);

fplot(y, [4.257475e-4 4.25749e-4]); ylim([-1.001 -1]); title('zoom fplot numeric')

fplot(eiY, [4.257475e-4 4.25749e-4]); ylim([-1.001 -1]); title('zoom fplot symbolic')

X = linspace(4.257475e-4, 4.25749e-4,1000);
YY = y(X);
subplot(2,1,1); plot(X, real(YY)); title('fixed numeric evaluation real')
subplot(2,1,2); plot(X, imag(YY)); title('fixed numeric evaluation imag')

X5idx = find(imag(YY), 5);
X5 = X(X5idx)

X5 = 1×5

0.000425747501501501 0.000425747504504504 0.000425747513513513 0.000425747516516516 0.00042574751951952

Y5 = YY(X5idx)

Y5 = 
          -1.00085203903026 +  1.38777878078145e-17i          -1.00085203903628 -  1.38777878078145e-17i          -1.00085203905432 -  1.38777878078145e-17i          -1.00085203906033 +  1.38777878078145e-17i          -1.00085203906635 -  1.38777878078145e-17i

syms t
subs(eiY, t, X5(:))

ans = 

vpa(ans)

ans =