# Inhomogeneous 2nd-order linear differential equation

I need to solve this:

$y'' + ay' + by = x(t)$

where nothing about the form of $x$ is known, except that it is bounded and non-negative. In addition it is known that $y(0) = 0$ and $y'(0) = 0$ (and $y''(0) = 0$ as well). I plugged this into the Wolfram ODE solver and got back a 'possible Lagrangian', which I really don't know how to use. (I have an undergraduate degree in math but I never took differential equations.) I also found this resource pt1, pt2, but it assumes that the solution is exponential in nature, which seems to violate my initial conditions since exponentials are never 0. I'm hoping there's a way to get my head around this without putting myself through the full course of diff eq so any help would be much appreciated.

Update: So actually Wolfram Alpha will solve this equation (simplified slightly in this edit), it was just the 'widget' that stopped short. However, it, as well as the techniques given in the references, explicitly assumes that the solution has an exponential form. Apparently this is a standard assumption but it is not the case for my data. $y$ is roughly sigmoidal. Coding up Wolfram's solution and choosing arbitrary constants $c_1$, $c_2$ so as to minimize the error in a case where $y$ was known produced a clearly wrong answer.

So, I will sharpen my question: How to solve the above differential question without assuming that $y$ is proportional to $e^{\lambda t}$?

• what are the $a,b,c$? – janmarqz Feb 3 '14 at 3:47
• @janmarqz They are arbitrary, typically values in $(0, 4]$. The equation derives from a model for $x$ where everything on the LHS is known; I'm essentially trying to invert the model. – Matt Phillips Feb 3 '14 at 3:53
• Just because there are exponentials in the answer doesn't mean you can't get $y(0)=0$. For example $y=e^t -e^{-t}$. The stuff in those resources is fine, but you can safely skip over the parts where there aren't constant coefficients. Maybe you'd be better off looking specifically for a resource on constant coefficient 2nd order linear ODEs. – mathematician Feb 3 '14 at 4:11
• Laplace transforms. – Chinny84 Feb 5 '14 at 16:16

I use the Laplace transform (see e.g. here ). Taking the Laplace transform of both sides we obtain $$(as^2+bs+c)Y(s)=X(s),$$ $Y(s),X(s)$ are the Laplace transform of $y(t),x(t)$, respectively. (In fact, $x(t)$ can be distribution also, the so-called generalized function.) So we have $$Y(s)=\frac{X(s)}{as^2+bs+c}.$$ Using the convolution theorem we obtain $$y(t)=x(t)*\mathcal{L}^{-1}\left(\frac{1}{as^2+bs+c} \right),$$ where $(f*g)(t)=\int_0^t f(t-v)g(v)\,dv$, and $f,g$ are piecewise continuous on $[0,\infty)$ (but this assumption can be relaxed).
Using the partial fraction decomposition $$\frac{1}{as^2+bs+c}=\frac{1}{a}\left(\frac{c_1}{s-s_1}+\frac{c_2}{s-s_2} \right),$$ or $$\frac{1}{as^2+bs+c}=\frac{1}{a(s-s_1)^2}.$$ Here $$\mathcal{L}^{-1}\left(\frac{1}{s-A} \right)=e^{At}, \qquad\mathcal{L}^{-1}\left(\frac{1}{(s-A)^2} \right)=e^{At}t.$$

• This looks like a solid way to go but I didn't understand $\frac{1}{as^2+bs+c}=\frac{1}{a(s-s_1)^2}$--you have three parameters, $a$, $b$, $c$ on the left and two $a$, $s_1$ on the left, how does that work? The use of Laplace transforms seems really promising though. – Matt Phillips Feb 5 '14 at 20:41
• @MattPhillips Please pay attention to "or". A quadratic equation has two roots. If they are differ ($s_1\neq s_2$), then we obtain the first decomposition, if they are equal ($s_1=s_2$) then we obtain the second one. – vesszabo Feb 6 '14 at 7:59

Let $f_\pm (t) := \exp(\lambda_\pm t)$ where $\lambda_\pm := (-a \pm \sqrt{a^2 - 4 b})/2$ be the fundamental solutions to the homogenous equation. Then the solution to the ODE, subject to $y(0)=r_0$ and $y^{'}(0)= \mu_0$, reads: \begin{equation} y(t) = \left(\begin{array}{cc} f_+(t) & f_-(t) \end{array}\right) \cdot \left( \begin{array}{cc} f_+(0) & f_-(0) \\ f_+^{'}(0) & f_-^{'}(0) \end{array} \right)^{-1} \cdot \left( \begin{array}{c} r_0 \\ \mu_0 \end{array} \right) + \int\limits_0^t \left(\frac{f_+(\xi) f_-(t)-f_-(\xi) f_+(t)}{f_+(\xi) f_-^{'}(\xi) - f_-(\xi) f_+^{'}(\xi)}\right) x(\xi) d\xi \end{equation} This can be proven using the Green's function method, for instance. Now, the first term on the right hand side is the generic solution to the homogenous equation and the second term is a specific solution to the inhomogenous equation.The quantity in the denominator in the integrand is called the Wronskian $W(\xi)$ and it reads: \begin{equation} W(\xi) = (\lambda_- - \lambda_+) e^{(\lambda_- + \lambda_+) \xi} = - \sqrt{a^2 - 4 b} e^{-a \xi} \end{equation} In our case $r_0=0$ and $\mu_0 = 0$ and so the solution reads: \begin{equation} y(t) = \frac{-1}{\sqrt{a^2 - 4 b}} \left( e^{\lambda_- t} \int\limits_0^t e^{(\lambda_++a) \xi} x(\xi) d \xi - e^{\lambda_+ t} \int\limits_0^t e^{(\lambda_-+a) \xi} x(\xi) d \xi \right) \end{equation}

• Great, I will check this out, although it looks very similar to what Wolfram produced, particularly in building the solution out of exponentials, but perhaps I just misunderstand the significance of assuming $f := exp(\lambda t)$. It's particularly nice that your solution contains no arbitrary constants, so +1 for an improvement on Wolfram Alpha, at least. – Matt Phillips Feb 5 '14 at 16:58
• The solution the I have quoted does not involve any arbitray constants because those constants have already been eliminated by using the initial conditions $y(0)=r_0$ and $y'(0)=\mu_0$. Another thing is that this solution satisfies any second order linear Ordinary Differential Equation(ODE), not only the one that you have quoted. On the other hand the Laplace transform method, despite its elegance usualy does not work if the coefficients of the ODE are not constant in time. – Przemo Feb 6 '14 at 10:54
• But what is the solution if $\;a^2 - 4 b = 0$ ? – Han de Bruijn Feb 11 '14 at 18:30

It's simple with Operator Calculus . Quite analogous to Example 1 on page 3 of this PDF we have: $$y'' + a y' + b y = x(t)$$ Operator Calculus enables us to "abstract" as much as possible from the solution $y(t)$ : $$\left[ \left(\frac{d}{dt}\right)^2 + a \frac{d}{dt} + b \right] y(t) = x(t)$$ What we are going to do next is decompose into factors. We are (almost) forced to do so, because the operator $(d/dt)$ and the constants $a$ and $b$ mutually behave as if they were ordinary numbers: the commutator of a differentiation and a constant is zero. But the commutative law for ( linear ) operators is the only thing which could be deviant from algebra with ordinary numbers. $$\left(\frac{d}{dt} \right)^2 + a \left(\frac{d}{dt} \right) + b = \left(\frac{d}{dt} - \lambda_1 \right) \left(\frac{d}{dt} - \lambda_2 \right)$$ $\lambda_1$ and $\lambda_2$ are roots of the so-called characteristic equation : $$\lambda^2 + a \lambda + b = 0 \quad \Longrightarrow \quad \lambda_{12}= (-a \pm \sqrt{a^2 - 4 b})/2$$ Herewith, the differential equation can be rewritten as: $$\left(\frac{d}{dt} - \lambda_1 \right) \left(\frac{d}{dt} - \lambda_2 \right) y(t) = x(t)$$ We are going to use now the "very useful formula" from Operator Calculus (please find it on page 2 of the abovementioned document) : $$\frac{d}{dt} + f = e^{-\int f\,dt}\; \frac{d}{dt} \; e^{+\int f\,dt}$$ In our case: $$\frac{d}{dt} - \lambda_{1,2} = e^{ - \int - \lambda_{1,2} \, dt} \ \frac{d}{dt}\ e^{ + \int - \lambda_{1,2} \, dt} = e^{\, \lambda_{1,2} t } \ \frac{d}{dt}\ e^{\, - \lambda_{1,2} t }$$ Giving, at last, for the O.D.E. : $$e^{ \lambda_1 t } \frac{d}{dt}\ e^{ - \lambda_1 t }\; e^{ \lambda_2 t } \frac{d}{dt}\ e^{ - \lambda_2 t } \; y(t) = x(t)$$ Systematic integration is possible now: $$e^{-\lambda_1 t} e^{\lambda_2 t} \frac{d}{dt}\ e^{-\lambda_2 t} y(t) = \int x(t) e^{ - \lambda_1 t }\,dt$$ $$y(t) = e^{\lambda_2 t} \int e^{(\lambda_1-\lambda_2) t} \left[ \int x(t) e^{ - \lambda_1 t }\,dt \right] \, dt$$ As has been said, $\lambda$ can be solved from $\lambda^2 + a\lambda + b = 0$, a quadratic equation with discriminant: $D=a^2-4b$. Real valued as well as complex solutions are possible. A special case for $D=0$ may be distinguished; in this case $\lambda_1 = \lambda_2$ and $e^{(\lambda_1-\lambda_2)t} = 1$ .
EDIT. Almost forgot the boundary conditions $y(0)=y'(0)=0$ to take into account: $$y(t) = e^{\lambda_2 t} \int_0^t e^{(\lambda_1-\lambda_2) u} \left[ \int_0^u x(v) e^{ - \lambda_1 v }\,dv \right] \, du$$

Note. If "nothing about the form of x is known" - but why a bounded & non-negative requirement - then $y''(0) = 0$ is in general not true , unless $x(0) = 0$ : $$y''(t) = \lambda_2^2 e^{\lambda_2 t} \int_0^t e^{(\lambda_1-\lambda_2) u} \left[ \int_0^u x(v) e^{ - \lambda_1 v }\,dv \right] \, du + (\lambda_1 + \lambda_2) e^{\lambda_1 t} \int_0^t x(v) e^{ - \lambda_1 v }\,dv \; + \; x(t)$$ With differential equations of the kind, $\,y''(0) = 0\,$ cannot serve as a boundary condition anyway.

How about using variation of parameters?

If you name $y_i(t)$ as the two parts of the homogeneous solution, obtained by assuming exponentials, so $$y(t) = A y_1(t) + B y_2(t) + y_p(t),$$ where $y_p(t)$ is the particular solution of the ODE, then by using variation of parameters as follows:

$$y(t) = A(t) y_1(t),$$

you will conclude that:

$$y_p(t) = y_1(t) \int M(t) \, dt,$$

where $M(t)$ is given by the following expression:

$$M(t) = \frac{e^{-at}}{y_1^2(t)} \int e^{at} y_1(t) x(t) \, dt.$$

You should apply now the initial conditions to the final solution, $y(t)$.

If you want me to elaborate some more on my answer, please let me know.

Cheers!