# Determining initial conditions for 2nd order non-homogenous differential equation

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}$$

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:

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:

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?

• 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. – Dylan Jun 28 '18 at 8:28
• 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. – Dylan Jun 28 '18 at 8:30
• @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. – a concerned citizen Jun 28 '18 at 9:01
• Is the ODE in $y$? If so, what is $x$? That's the part that confused me. – Dylan Jun 28 '18 at 9:39
• @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. – a concerned citizen Jun 28 '18 at 9:47

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

• 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)? – a concerned citizen Jun 28 '18 at 11:22
• 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? – a concerned citizen Jun 28 '18 at 11:27
• 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. – Lutz Lehmann Jun 28 '18 at 11:28
• 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}$. – Lutz Lehmann Jun 28 '18 at 11:32

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}

• 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... :-) – a concerned citizen Jun 28 '18 at 12:40