5 views (last 30 days)

Show older comments

Walter Roberson
on 3 Aug 2021

You do not seem to be able to find a closed form solution with symbolic d or e.

Your ilaplace() involves the sum of an expression involving roots of a cubic function. Two of the roots are (probably) complex valued.

When I do some tests with the complex values, it looks as if the imaginary parts vanish. However, at the moment I do not see any way of proving that, as the coefficients of the imaginary terms are not small. I suspect you would have to prove that the phases cancel out because the two complex roots would have to be complex conjugates.

The symsum() of roots of a polynomial is a nuisance: mathworks does not define any way to expand the roots, just using root indices instead. Because it is only a cubic, it is possible to replace the indexed root with its closed form formula... getting a long formula, but something that has more potential to be integrated.

Substituting in closed form polynomial roots instead of the placeholders that ilaplace uses, is not trivial. This is the first time I have been able to put suitable code together, and I am not certain that I got it correct.

Taking sign(7.4513-z) implies that you think z is real-valued, which seems likely to be true in theory, but is difficult to prove, and might not be true numerically.

The ilaplace result would normally have a couple of derivatives of dirac deltas on t, and those derivatives of dirac(0) would normally contribute infinite values... it is not obvious to me that the ilaplace result has a well-defined value at t = 0. I skip the problem by adding an assumption that t is strictly greater than 0; in theory that is a problem for the integration starting from 0. The integration should probably start from realmin() instead of 0.

If you substitute in specific numeric values for d and e in the below, the int(w) might take a quite long time. You probably should not do that, and should probably use vpaintegral() instead of int() for that case.

Q = @(v) sym(v);

syms s t d e real

assume(t>0);

y = (Q(0.005659)*s^2+Q(0.05206)*s+Q(0.004496)/(s^3+Q(9.357)*s^2+Q(2.228)*s+Q(0.06147)))*(Q(27.5084)*d+Q(27.5084)*(e/s));

z = ilaplace(y,s,t);

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);

ez = simplify(expand_symsum(z));

Bound = Q(7.4513);

w = piecewise(ez < Bound, Bound-ez, ez-Bound);

I = int(w, t, 0, 2000);

w20 = vpa(simplify(w),20);

w20z = vpa(real(mapSymType(w20, 'vpareal', @(X7) piecewise(abs(X7)<1e-10,0,X7))),20)

I20 = int(w20z, t, 0, 2000)

Walter Roberson
on 14 Aug 2021

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

Start Hunting!