1
$\begingroup$

With the help of wxMaxima (math not being my language), I am trying to find out the impulse response for a 2nd order transfer function:

$$H(s)=\frac{a_2s^2+a_1s+a_0}{b_2s^2+b_1s+b_0}=\frac{2.3s^2+0.1s+4.5}{1.3s^2+0.7s+1.1}$$

If I use the inverse Laplace, ilt(partfrac(H(s), s), s, t);, it gives the correct response:

$$\exp\left(-\frac{7 t}{26}\right)\left( \frac{9668 \sin{\left( \frac{\sqrt{523} t}{26}\right) }}{169 \sqrt{523}}-\frac{148 \cos{\left( \frac{\sqrt{523} t}{26}\right) }}{169}\right)$$

(I omitted the ilt(a2/b2,s,t) term since it cannot be plotted) but if I am trying to implement H(s) into a differential equation:

$$b_2y''(t)+b_1y'(t)+b_0y(t)=a_2x''(t)+a_1x'(t)+a_0x(t)$$ $$1.3y''(t)+0.7y'(t)+1.1y(t)=2.3x''(t)+0.1x'(t)+4.5x(t)$$

and use wxMaxima's builtin ode2() and ic2() to solve it:

b2*'diff(y, t, 2) + b1*'diff(y, t) + b0*y = a2*'diff(x, t, 2) + a1*'diff(x, t) + a0*x;
ode2( ''%, y, t);
ic2(%th(1), t=0, y=0, 'diff(y,t)=1);

the result is different. I used for the initial conditions $y'$=1, since the input is the Dirac function and I cannot make it $\infty$, and $y$=0, and the result is this:

$$\exp\left(-\frac{7 t}{26}\right)\left( -\frac{\sin{\left( \frac{\sqrt{523} t}{26}\right) } \left( 315 \sqrt{523} x-279 \sqrt{523}\right) }{5753}-\frac{\cos{\left( \frac{\sqrt{523} t}{26}\right) } \left( 45 x+1\right) }{11}\right) +\frac{45 x+1}{11}$$

r1

The red trace is the inverce Laplace. Both traces are missing the initial Dirac because it can't be plotted. For the same reason, the ODE is plotted with x=0, and also since any sort of impulse (e.g. if t<0.1 then 1 else 0) or similar don't work.

If I use an initial output value for $y$ of -148/169 (which is the cos() term of the inverse Laplace) and I further tweak the $y'$ initial condition to be 2.33 (with the hammer), I get this:

r2

It looks like there is a DC component that just stays there, no matter what. Since my intuition fails, maybe the hammer won't, so I thought of differentiating the answer with changed initial conditions, which would give zero steady-state. Here's with $y$=-2.22 and $y'$=-148/169:

r3

but this is some brute-force that I don't think it applies.

My questions: Is what I am doing the correct way of solving the ODE? If no, how? If yes, how should I determine the correct initial conditions?

$\endgroup$
  • $\begingroup$ Since $H(s\to \infty) \ne 0$, you'll need to take out the constant, i.e. $$ H(s) = \frac{a_2}{b_2} + \frac{c_1s + c_0}{b_2s^2+b_1s+b_0} $$ So there needs to be added term of $\frac{a_2}{b_2}\delta(t)$ in the inverse. $\endgroup$ – Dylan Jun 28 '18 at 8:28
  • $\begingroup$ I'm also not quite sure I understand your question. Are you looking for the initial conditions of an IVP where $H(s)$ is the Laplace transform of the solution? Or is it a different IVP where $H(s)$ is just the transform of the inhomogeneous term? Please clarify by writing down the actual equation. $\endgroup$ – Dylan Jun 28 '18 at 8:30
  • $\begingroup$ @Dylan partfrac() inside ilt() takes care of that, but the builtin delta(), or ilt(a2/b2,s,t) cannot be plotted, so I eliminated it. The equation is a bit obscured inside the wxMaxima lines, I'll add it in clear. I am confused as to how to solve the differential equation so when I plot both the results of the inverse Laplace and the result of the ODE, they should match. It seems to be related to the initial conditions, and yet not quite: the 1st plot (blue) shows the brute result, with $y$=0 and $y'$=1, the 2nd with tweaked IVP, 3rd with both tweaked and the derivative of the result. $\endgroup$ – a concerned citizen Jun 28 '18 at 9:01
  • $\begingroup$ Is the ODE in $y$? If so, what is $x$? That's the part that confused me. $\endgroup$ – Dylan Jun 28 '18 at 9:39
  • $\begingroup$ @Dylan Ah, I'm sorry, $y$ is the output, $x$ is the input, while $t$ is the variable. In this case, $x$ is the Dirac function, but I'm hoping that the result be true for any input. $\endgroup$ – a concerned citizen Jun 28 '18 at 9:47
1
$\begingroup$

The general setup is that all functions are zero for $t<0$. You get some difficulties as the right side contains the second derivative of the Dirac delta distribution if $x=\delta$. Let's reduce the derivatives on the right side by shifting $y$ by a corresponding multiple of $x$. Define $$u=y-\frac{a_2}{b_2}x,$$ then $$ b_2u''+b_1u'+b_0u=\frac{b_2a_1-b_1a_2}{b_2}x'+\frac{b_2a_0-b_0a_2}{b_2}x $$ which tells us that to compensate for the next highest order distributional term, $u$ has to contains a part that is an anti-derivative $X$ of $x$, $X'=x$. For $x=\delta$ it is a unit ramp. Let $$v=u-\frac{b_2a_1-b_1a_2}{b_2^2}X,$$ then $$ b_2v''+b_1v'+b_0v=\frac{b_2^2a_0-b_0b_2a_2-b_1b_2a_1+b_1^2a_2}{b_2^2}x-\frac{b_0b_2a_1-b_0b_1a_2}{b_2^2}X $$ To remove the last singular part, compensate with an anti-derivative $\Xi$ of $X$, $\Xi'=X$, $\Xi''=x$. Let $$w=v-\frac{b_2^2a_0-b_0b_2a_2-b_1b_2a_1+b_1^2a_2}{b_2^3}\Xi.$$ Then we know that $w$ satisfies an ODE with piecewise continuous right side, is thus a continuously differentiable function (piecewise $C^2$). From the general condition $w(t)=0$ for $t<0$ we get the initial conditions $w(0)=0$, $w'(0)=0$ and in combination of all substitutions using $x=\delta$, $X(t)=t_+^0=1$ for $t>0$, $\Xi(t)=t_+=\max(0,t)$ $$ y(t)=w(t)+\frac{a_2}{b_2}\delta(t)+\frac{b_2a_1-b_1a_2}{b_2^2}(t_+)^0+\frac{b_2^2a_0-b_0b_2a_2-b_1b_2a_1+b_1^2a_2}{b_2^3}t_+. $$ Apart from the $δ$ singularity around $t=0$ this has an initial value $$\lim_{t\to +0}y(t)=\frac{b_2a_1-b_1a_2}{b_2^2}$$ and $$\lim_{t\to +0}y'(t)=\frac{b_2^2a_0-b_0b_2a_2-b_1b_2a_1+b_1^2a_2}{b_2^3}.$$ Inserting the numerical values this gives $$ y(0_{+0})=-0.8757396449704141\\ y'(0_{+0})=2.4360491579426493 $$ which gives the solution graph enter image description here Python code

a2,a1,a0,b2,b1,b0  =  2.3, 0.1, 4.5, 1.3, 0.7, 1.1
y0 = (b2*a1-b1*a2)/b2**2
y1 = (b2**2*a0-b0*b2*a2-b1*b2*a1+b1**2*a2)/b2**3
z = (-b1+(b1**2+0j-4*b0*b2)**0.5)/(2*b2)
a, b = z.real, z.imag
t = np.linspace(0,25,1501)
m = (y1-a*y0)/b;# y1=y'(0)=a*y0+m*b 
y = np.exp(a*t)*(y0*np.cos(b*t)+m*np.sin(b*t)) 
plt.plot(t,y); plt.grid(); plt.show()
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the answer, it certainly brings a new way to look at it. Still, in wxMaxmia, I still get the extra $(45x+1)/11$ term, and this is what I plot with the new i.c.: i.stack.imgur.com/1HaVO.png (blue=ODE). If I remove the offending term: i.stack.imgur.com/Ky6Sk.png (with i.stack.imgur.com/GcZrw.png being the difference). Maybe wxMaxima has problem interpreting x as $\delta$(t)? $\endgroup$ – a concerned citizen Jun 28 '18 at 11:22
  • $\begingroup$ Also, for a step input, would I have to do the integration only once? That is, once I get the ramp, there's no need to go further, since the function is continuous? Or I have to keep going twice, so that the resulting equation is, again 2nd degree? $\endgroup$ – a concerned citizen Jun 28 '18 at 11:27
  • $\begingroup$ That is certainly possible, as for $t>0$ the equation is simply the homogeneous $b_2y''+b_1y'+b_0y=0$ which does not have non-trivial constant solutions. $\endgroup$ – Lutz Lehmann Jun 28 '18 at 11:28
  • $\begingroup$ Yes, for $x$ the unit step function $(t_+)^0$ you would get that $v$ is already $C^1$ and piecewise $C^2$ so that $v(0)=v'(0)=0$ and then $y(t+0)=\frac{a_2}{b_2}$ and $y'(t+0)=\frac{b_2a_1-b_1a_2}{b_2^2}$. $\endgroup$ – Lutz Lehmann Jun 28 '18 at 11:32
1
$\begingroup$

A more systematic way considers a solution $y_0$ of $$L(y_0)=b_2y_0''+b_1y_0'+b_0y_0=x,$$ $y(t)=0$ for $t<0$. Then \begin{align} y&=a_0y_0+a_1y_0'+a_2y_0''\\ &=\left(a_0-\frac{a_2b_0}{b_2}\right)y_0+\left(a_1-\frac{a_2b_1}{b_2}\right)y_0'+ \frac{a_2}{b_2}x \end{align} is a solution for the given ODE with derivative \begin{align} y'&=\left(a_0-\frac{a_2b_0}{b_2}\right)y_0'+\left(a_1-\frac{a_2b_1}{b_2}\right)y_0''+ \frac{a_2}{b_2}x'\\ &=-\frac{b_0}{b_2}\left(a_1-\frac{a_2b_1}{b_2}\right)y_0+\left(a_0-\frac{a_2b_0+a_1b_1}{b_2}+\frac{a_2b_1^2}{b_2^2}\right)y_0'+\left(\frac{a_1}{b_2}-\frac{a_2b_1}{b_2^2}\right)x+\frac{a_2}{b_2}x' \end{align}


For $x=\delta$ we get $y_0=uv$ where $u(t)=0$ for $t<0$, $u(t)=1$ for $t>0$ is the unit jump and $v$ a solution to $L(v)=0$ with $v(0)=0$, $v'(0)=\frac1{b_2}$, which implies $y_0'=δv+uv'=δv(0)+uv'=uv'$. Thus \begin{alignat}{1} y(0_{+0})&=\left(a_1-\frac{a_2b_1}{b_2}\right)v'(0)&=\frac{a_1}{b_2}-\frac{a_2b_1}{b_2^2} \\ y(0_{+0})&=\left(a_0-\frac{a_2b_0+a_1b_1}{b_2}+\frac{a_2b_1^2}{b_2^2}\right)v'(0)&=\frac{a_0}{b_2}-\frac{a_2b_0+a_1b_1}{b_2^2}+\frac{a_2b_1^2}{b_2^3} \end{alignat}

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the extra details. It's a bit hard to chew for my teeth, but I'll get it done, however, I have to say that this topic seems to be something you really don't like talking about, at all... :-) $\endgroup$ – a concerned citizen Jun 28 '18 at 12:40

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.