5
$\begingroup$

I am having some trouble when using OutputResponse. For instance consider the following command

Simplify[Chop[OutputResponse[
  TransferFunctionModel[
      {
        {{12265.875860667435 + 15863.964729950849*s + 5000.*s^2}}, 
        12265.875860667435 + 15986.623488557521*s + 5132.222766045224*s^2 +
        81.0807167110999*s^3 + 16.908502769264427*s^4 + 1.*s^5
      },
      s,
      SystemsModelLabels -> {{None}, {None}}
      ],
  UnitStep[t],
  t]]]

The output is zero on version 9.0.1-mac. It can't be zero.

The above transfer function is a result of the feedback connection of a series connection between two transfer functions.

What am I missing?

Many thanks

Ed

Improve this question
$\endgroup$
2
  • $\begingroup$ The unChop'd result has products comprised of very small times very large factors. Use Expand before you use Chop so that these will collapse to numbers of reasonable size. And hope you don't get bitten hard by cancellation or truncation error. $\endgroup$
    – Daniel Lichtblau
    Feb 4, 2013 at 15:15
  • $\begingroup$ @Daniel - Many thanks. That seems to work. I am puzzled at how to use the functions that comes with the Control Package. Something as simple as closing the loop and then calculating the output response can lead to the problems I am experiencing. $\endgroup$
    – Ed Mendes
    Feb 5, 2013 at 17:53

3 Answers 3

Reset to default
7
$\begingroup$

For machine-precision computations, lots of coefficients $\sim 10^{-15}$ may be a sign that machine precision is actually not sufficient (and, yes, such results can be different in different versions or even on different platforms). The standard trick in such cases is to Rationalize[] the input or use a higher precision, at the expense of slower, but more reliable computations.

For instance, this is reproducible in V8 and V9 on my machine:

In[48]:= OutputResponse[
   Rationalize[
    TransferFunctionModel[{{{12265.875860667435 + 
         15863.964729950849*s + 5000.*s^2}}, 
      12265.875860667435 + 15986.623488557521*s + 
       5132.222766045224*s^2 + 81.0807167110999*s^3 + 
       16.908502769264427*s^4 + 1.*s^5}, s, 
     SystemsModelLabels -> {{None}, {None}}], 0], 1, t] // N // Chop

Out[48]= {0. + 
  6.3706*10^-7 (0. - 
     784856. (0.00285107 E^(-22.2373 t) - 4.59391 E^(-1.88023 t) + 
        6.60436 E^(-1.30805 t) - (0.00664625 - 
           0.00421948 I) E^((4.25855 - 14.3576 I) t) - (0.00664625 + 
           0.00421948 I) E^((4.25855 + 14.3576 I) t)))}

(I replaced UnitStep[t] -> 1, which is equivalent in this case, but avoids a lot of unnecessary symbolic work and runs much faster.)

Or, in V9, you can try:

In[50]:= OutputResponse[
   SetPrecision[
    TransferFunctionModel[{{{12265.875860667435 + 
         15863.964729950849*s + 5000.*s^2}}, 
      12265.875860667435 + 15986.623488557521*s + 
       5132.222766045224*s^2 + 81.0807167110999*s^3 + 
       16.908502769264427*s^4 + 1.*s^5}, s, 
     SystemsModelLabels -> {{None}, {None}}], 20], 1, t] // 
  Simplify // Chop

Out[50]= {0.2452185003 - 0.2477786589 E^(-22.2373307055086 t) + 
  0.02140659510 E^(-1.8802349880252 t) - 
  0.01884643653 E^(-1.3080450008032 t) - 
  0.7547814997 E^(4.2585539625363 t) Cos[14.3575707114421 t] + 
  0.7547814993 Cos[14.3575707114421 t]^2 - 
  0.1588054718 E^(4.2585539625363 t) Sin[14.3575707114421 t] + 
  0.7547814997 Sin[14.3575707114421 t]^2}

And finally, when you only need to plot an output response, or do something similar, you may want to use the numeric syntax OutputResponse[sys, u, {t, 0, tmax}] rather than the symbolic one, OutputResponse[sys, u, t]. The numeric route (in this case, via NDSolve) can be faster, avoids many pitfalls, and is also reproducible between V8 and V9:

In[54]:= OutputResponse[
 TransferFunctionModel[{{{12265.875860667435 + 15863.964729950849*s + 
      5000.*s^2}}, 
   12265.875860667435 + 15986.623488557521*s + 
    5132.222766045224*s^2 + 81.0807167110999*s^3 + 
    16.908502769264427*s^4 + 1.*s^5}, s, 
  SystemsModelLabels -> {{None}, {None}}], UnitStep[t], {t, 0, 1}]

Out[54]= {12265.875860667435*InterpolatingFunction[][t] + 
  15863.964729950849*InterpolatingFunction[][t] + 
     5000.*InterpolatingFunction[][t]}

In[55]:= Plot[%, {t, 0, 1}]
...
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for pointing out various options. The help pages on Mathematica documentation do not seem to cover any of them. $\endgroup$
    – Ed Mendes
    Feb 19, 2013 at 22:51
1
$\begingroup$

I don't know what the output should be, but it's not zero if you remove the Chop:

Simplify[OutputResponse[TransferFunctionModel[{{{12265.875860667435 + 15863.964729950849*s + 5000.*s^2}}, 12265.875860667435 + 15986.623488557521*s + 5132.222766045224*s^2 + 81.0807167110999*s^3 + 16.908502769264427*s^4 + 1.*s^5}, s,SystemsModelLabels -> {{None}, {None}}], UnitStep[t], t]]

{E^(-55.1098 t) ((8.52348*10^-17 - 7.19335*10^-32 I) E^(
     32.8724 t) + (7.98901*10^-18 + 1.47911*10^-31 I) E^(
     53.2295 t) - (8.84652*10^-18 - 1.47911*10^-31 I) E^(
     53.8017 t) + (4.2727*10^-16 - 1.64474*10^-31 I) E^(59.3683 t)
      Cos[14.3576 t] + (1.83851*10^-16 + 4.76668*10^-32 I) E^(
     59.3683 t)
      Sin[14.3576 t] + ((-0.247779 + 8.38088*10^-17 I) E^(
        32.8724 t) + (4.10452*10^-16 - 4.40027*10^-31 I) E^(
        34.1805 t) - (1.35081*10^-15 - 1.05004*10^-17 I) E^(
        34.7527 t) + (0.0214066 + 4.44089*10^-16 I) E^(
        53.2295 t) - (0.0188464 + 6.66134*10^-16 I) E^(
        53.8017 t) - (1.8673*10^-15 + 3.15544*10^-30 I) E^(
        54.5376 t) + (0.245219 + 8.88178*10^-16 I) E^(
        55.1098 t) - (4.12448*10^-16 + 2.67679*10^-17 I) E^(
        55.682 t) - (4.30823*10^-17 - 1.55326*10^-31 I) E^(
        75.4669 t) + (1.03401*10^-17 + 5.54668*10^-32 I) E^(
        76.0391 t) + ((-3.33118*10^-17 - 1.0091*10^-16 I) E^(
           28.6139 t) - (1.64219*10^-17 - 8.98175*10^-18 I) E^(
           48.971 t) - (1.79284*10^-17 - 1.15633*10^-18 I) E^(
           49.5432 t) - (0.754781 - 3.99637*10^-16 I) E^(
           59.3683 t) + (1.40072*10^-15 + 8.95624*10^-17 I) E^(
           60.6764 t) - (4.62884*10^-15 + 1.75057*10^-16 I) E^(
           61.2486 t) + (1.13964*10^-17 + 7.12682*10^-18 I) E^(
           81.6057 t)) Cos[
         14.3576 t] + (0.754781 - 4.01372*10^-16 I) E^(55.1098 t)
         Cos[14.3576 t]^2 + ((3.58981*10^-16 + 1.69519*10^-17 I) E^(
           28.6139 t) + (8.3687*10^-18 + 3.06739*10^-18 I) E^(
           48.971 t) + (1.06871*10^-17 + 8.15728*10^-19 I) E^(
           49.5432 t) - (0.158805 - 8.34836*10^-17 I) E^(
           59.3683 t) + (5.50678*10^-16 + 2.29738*10^-16 I) E^(
           60.6764 t) - (1.97153*10^-15 + 3.1209*10^-16 I) E^(
           61.2486 t) + (1.4127*10^-16 - 1.3553*10^-17 I) E^(
           81.6057 t)) Sin[
         14.3576 t] + (0.754781 - 6.16044*10^-16 I) E^(55.1098 t)
         Sin[14.3576 t]^2 + (4.55706*10^-16 + 1.7103*10^-16 I) E^(
        55.1098 t) Sin[28.7151 t]) UnitStep[t])}

Another option is to perform the Simplify first and then do the Chop, which yields the following:

{E^(-55.1098 t) (-0.247779 E^(32.8724 t) + 0.0214066 E^(53.2295 t) - 
    0.0188464 E^(53.8017 t) + 0.245219 E^(55.1098 t) - 
    0.754781 E^(59.3683 t) Cos[14.3576 t] + 
    0.754781 E^(55.1098 t) Cos[14.3576 t]^2 - 
    0.158805 E^(59.3683 t) Sin[14.3576 t] + 
    0.754781 E^(55.1098 t) Sin[14.3576 t]^2) UnitStep[t]}

Plot[%,{t,0,1}]

enter image description here

Improve this answer
$\endgroup$
5
  • $\begingroup$ Many thanks. Out of curiosity what version are you running? The reason for the question comes to the fact that I have got the same results as the ones you have shown in your reply when I run Mathematica Version 8.0.4 and issue the exact same command as the one in my question. $\endgroup$
    – Ed Mendes
    Feb 3, 2013 at 9:25
  • $\begingroup$ I'm also running v9.0.1 on Mac OS X and I get the same results you describe when I run your example on MMA v8. Without spending too much time on this, it looks to me like MMA produces equivalent but different symbolic results under 8 and 9, and the differences become important when you apply Chop. For that reason (and many others), you should exercise caution when "Chopping". I generally only use Chop when I'm trying to display numbers in a "clean" manner, but never in the midst of a calculation. $\endgroup$
    – Cassini
    Feb 3, 2013 at 13:09
  • $\begingroup$ Funny the results in MMA v8 are definitely different from those of MMA v9 (even if I use or don't use CHOP). It is getting more and more difficult to accept the results that come from the Control Package. Forget mention that if I do not use CHOP on v9 the output is huge with lots and lots of 10ˆ(-15) $\endgroup$
    – Ed Mendes
    Feb 3, 2013 at 16:35
  • $\begingroup$ @EdMendes: FYI, I've noticed you're fairly new to this forum and that you haven't "accepted" any answers to your questions. If you receive useful answers to your questions, you should vote them up. If you receive an answer which truly satisfies your question, you should "accept" it. $\endgroup$
    – Cassini
    Feb 3, 2013 at 16:50
  • $\begingroup$ Many thanks for pointing that out. $\endgroup$
    – Ed Mendes
    Feb 3, 2013 at 19:10
1
$\begingroup$

A few observations:

  1. Perhaps I'm missing something, but after reviewing the output one can see that all the terms are of the order 10^-15 or smaller. I believe that using the command Chop will reduce all of these to the integer zero.
  2. The impulse response (with DiracDelta input) appears to show an unstable system in v 9.0.1 (see results below).

Using version 9.0.1 I get the following results, which differ from those already posted:

a. The response to a UnitStep input:

responseUnit = 
Simplify[OutputResponse[
TransferFunctionModel[{{{12265.875860667435 + 15863.964729950849*s + 5000.*s^2}}, 12265.875860667435 + 15986.623488557521*s + 5132.222766045224*s^2 + 81.0807167110999*s^3 + 16.908502769264427*s^4 + 1.*s^5}, s, 
SystemsModelLabels -> {{None}, {None}}],
UnitStep[t], t]];

gives me the following output:

OutputResponse to UnitStep input for 0 < t < 1

b. The response to a DiracDelta input:

responseUnit = 
Simplify[OutputResponse[
TransferFunctionModel[{{{12265.875860667435 + 15863.964729950849*s + 5000.*s^2}}, 12265.875860667435 + 15986.623488557521*s + 5132.222766045224*s^2 + 81.0807167110999*s^3 + 16.908502769264427*s^4 + 1.*s^5}, s, 
SystemsModelLabels -> {{None}, {None}}],
DiracDelta[t], t]];

gives me the following output: enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Many thanks. Notice that b) contradicts zero as the response of the system to an unit step. Regarding a) Have you noticed that the response vanishes after 0.3 secs. That does not suppose to happen. And yes, the system is unstable but that does not mean one can not calculate the output response. $\endgroup$
    – Ed Mendes
    Feb 3, 2013 at 16:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.