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

Question

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?

Answer 1

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

Answer 2

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}