Higher-order linear differential equations

In the previous topic, we learned about an important special case of ODEs: first-order linear differential equations. Despite their prevalence, many interesting problems in various fields of science are described by higher-order differential equations. We are now going to turn our focus to them.

Definition

Let's recall the concept of linearity from the previous topic. All the derivatives $y^{(n)}(x)$ will enter our equation in the first degree.

In order to keep things simple, in this topic, we will consider all the coefficients of the derivatives to be constants. We would be happy to consider a general equation of the form

$$a_n(x)y^{(n)}(x) + a_{n-1}(x)y^{(n-1)}(x) + ... + a_1(x)y'(x) =f(x)$$...but, unfortunately, it is impossible to solve such an equation in the general form. This is a common situation in the field of differential equations: as soon as you take a step off the road paved by theory, you lose the opportunity to solve the equation analytically. Fortunately, there are numerous methods for solving ODEs numerically, but we'll talk about them some other time.

The equations we will deal with today will have the form

$$a_ny^{(n)}(x) + a_{n-1}y^{(n-1)}(x) + ... + a_1y'(x) =f(x), \\ a_n, a_{n-1},..., a_1 = const$$Equations of degree greater than 2 are rarely encountered in practice. In order not to overload the text with indices, we will consider a particular case and then see how it can be extended to equations of higher orders.

To sum up, we will devote most of our presentation to equations of this kind:

$$y''(x)+py'(x) + qy(x) = f(x),\\ p,q = const \tag{*}$$Recall that an equation with $f(x) \equiv 0$ is called homogeneous. Obviously, every nonhomogeneous equation has a homogeneous counterpart. Solving this counterpart is a crucial step in finding a solution for the initial problem.

Solving higher-order linear homogeneous ODEs

Let's consider a homogeneous equation:

$$y''(x)+py'(x) + qy(x) = 0,\\ p,q = const \tag{**}$$To solve it, we should construct the characteristic equation:

$$\lambda^2+p\lambda+q = 0$$The principle of its construction is obvious – instead of the derivative $y^{(n)}(x)$, we write lambda raised to the appropriate power: $\lambda^n$.

Why should we do that and what's the point of it? It turns out that the roots of the characteristic equation completely determine the solutions of the original differential equation $(**)$. We will not strictly prove this fact (the proof can be found in textbooks), but you can always substitute the found solutions and check whether they satisfy our homogeneous ODE.

Let's define the roots (no matter real or complex) of our characteristic equation as $\lambda_1, \lambda_2$.

As it always happens when solving quadratic equations, we have three cases:

1) Two different real roots (discriminant $D>0$)

In this case, the solution of the equation $(**)$ takes the form

$$y(x) = C_1e^{\lambda_1x} + C_2e^{\lambda_2x} \tag{1}$$

Recall that the solution to an $n$-th order ODE has $n$ arbitrary constants (in our case, $n=2$)!

2) Two different complex roots (discriminant $D<0$)

If the roots of a quadratic equation are complex, they must be complex conjugate roots. Let's present the roots in the form

$$\lambda_1 = \alpha + i\beta, \\ \lambda_2 = \alpha - i\beta, \\ \alpha, \beta \in \mathbb{R}$$The solution of $(**)$ in this case is:

$$y(x) = e^\alpha(C_1cos\beta x + C_2sin\beta x) \tag{2}$$Those of you who are familiar with complex exponents will easily understand that this case is virtually no different from the first one. We just have to substitute the two building blocks in the solution $(1)$ with their linear combinations:

$$y_1(x) = e^{\lambda_1x} = e^{(\alpha + i\beta)x} = e^\alpha (cos \beta x + isin\beta x)\\ y_2(x) = e^{\lambda_2x} = e^{(\alpha - i\beta)x} = e^\alpha (cos \beta x - isin\beta x)\\ y_1^* = \frac{y_1+y_2}{2} = e^\alpha cos \beta x\\ y_2^* = \frac{y_1-y_2}{2i} = e^\alpha sin \beta x\\ y(x) = C_1y_1^*(x) + C_2y_2^*(x) =e^\alpha(C_1cos\beta x + C_2sin\beta x)$$

3) Two repeated real roots (discriminant $D=0$)

This case is a little bit different. The solution of equation $(**)$ will be:

$$y(x) = e^{\lambda x}(C_1 + C_2x) \tag{3}$$where $\lambda = \lambda_1 = \lambda_2$.

A small remark about equations of degree greater than 2. As we know, any algebraic equation of degree $n$ has exactly $n$ complex roots (some or all of them may be real). Using these roots, we can construct a solution to any homogeneous equation just as we did for the case of $n=2$. Each unique real root gives an exponent in the solution (see point 1). Each pair of complex conjugate roots gives an exponent multiplied by the "trigonometric block" (see point 2). Each root repeated $k$ times gives a polynomial of degree $k-1$ multiplied by the exponent with that root (see point 3).

Numerical example

Let's see how this is done in practice. We will consider the following equation:

$$y'' -2y' + 10y = 0$$Its characteristic equation is written as follows:

$$\lambda^2-2\lambda+10 = 0\\$$It has two complex conjugate roots (case 2):

$$D = 2^2-4\cdot10 = -36 = (6i)^2\\ \lambda_{1,2} = \frac{-2 \pm 6i}{2} = -1 \pm 3i$$

Therefore, according to $(2)$, the solution of the initial equation is:

$$y(x) = e^{-x}(C_1cos3x + C_2sin3x), \\ C_1, C_2 = const$$

Higher-order linear nonhomogeneous ODEs

Now, let's turn our attention to nonhomogeneous equations. As before, the first step in solving such an equation will be solving its homogeneous counterpart. The whole algorithm is as follows:

1. Find a solution of a homogeneous equation $(**)$ $y_0(x)$.
2. Find some particular solution of the initial equation $(*)$ $y_p(x)$.
3. Construct a full solution of (*) by summation: $y(x) = y_0(x) + y_p(x)$. The result is called the general solution of $(*)$.

Step 1 has been discussed in detail above, and step 3 is self-explanatory. So now we will focus on step 2 – finding a particular solution.

We will use the method of undetermined coefficients to find a particular solution $y_p(x)$. First, we will "guess" the form of the solution, leaving the numerical coefficients in this form undefined. Then, we will substitute our ansatz into the differential equation $(*)$ and determine the values of the coefficients.

A legitimate question arises – how to guess the form of the particular solution? Unfortunately, it is impossible to give a single recipe here – it all depends on the type of function $f(x)$ in the right-hand side of the equation (*). However, there are many useful heuristics that allow you to "guess" the right form of the solution. Now we will discuss some of them.

• $f(x)$ is a polynomial, for example $f(x) = x^2 + 1$. Then $y_p(x)$ should be searched among polynomials of the same order multiplied by an additional term: $y_p(x) = (Ax^2 + Bx + C)\cdot x^z$, where $z$ is the number of zero roots of the characteristic equation. Yeah, it sounds convoluted and it actually is. But don't worry, in most second-order equations $\lambda_{1,2} \neq 0$ and $z=0$. Note that the coefficients should not be missed even if they are absent in the $f(x)$.
• $f(x)$ is an exponent, for example $f(x) = e^{4x}$ . If the coefficient in the exponent doesn't coincide with the roots of the characteristic equation ($\lambda_{1,2} \neq 4$), we should search for $y_p(x)$ in the same form: $y_p(x) = Ae^{4x}$. Otherwise, the ansatz should be multiplied by x: $y_p(x) = Axe^{4x}$.
• $f(x)$ is a sin or a cos function, for example $f(x) = sin3x$. The particular solution takes the form $y_p(x) = Asin3x + Bcos3x$. Note that it always contains both sin and cos parts, no matter if they are present in $f(x)$.

These cases cover most of the situations you may encounter in practice. However, another important aspect remains – what if f(x) consists of several of the "building blocks" discussed above? Here are two simple rules to help you figure out the form of $y_p(x)$ in this situation:

1. If $f(x)$ contains an exponential, we should ignore it and write down the guessed form of the solution for the remaining part. Then tack the exponential back on without any coefficient before it.
   
2. When $f(x)$ is a product of a trigonometric function with a polynomial, we first write down the ansatz for the polynomial only and multiply that by the appropriate cosine. Then we add the same polynomial with different coefficients and multiply that by the appropriate sine.

Let's have a look at a concrete example. Let $f(x) = e^{4x}(x^2+1)sin3x$. Following the rules above, we forget about the exponent and start constructing the template for the polynomial: $$Ax^2 + Bx+C$$Since we have a sin part, we should now concatenate this template with the trigonometric functions:

$$(Ax^2 + Bx+C)sin3x +(Dx^2+Ex+F)cos3x$$Finally, we return the exponent and obtain the guess for $y_p(x)$:

$$y_p(x) = e^{4x}((Ax^2 + Bx+C)sin3x +(Dx^2+Ex+F)cos3x)$$A final remark: if $f(x) = f_1(x) + f_2(x)$, where $f_1(x)$ and $f_2(x)$ belong to different types discussed above, the particular solution should be divided into blocks and found step by step, ignoring other parts. After that, the final particular solution may be found as a simple sum of the blocks: $y_p(x) = y_{p1}(x) + y_{p2}(x)$.

A complete example

This journey through theory was quite exhausting. Let's consider the whole process of solving an equation $(*)$ from the beginning to end with a concrete example. We will be guided by the algorithm provided in the previous section.

$$y'' -2y' + 10y = xe^{2x}$$

1. Step 1: the homogeneous equation. We've already solved it before:$$y_0(x) = e^{-x}(C_1cos3x + C_2sin3x), \\ C_1, C_2 = const$$
2. Step 2: finding a particular solution. The right-hand side contains a polynomial and the exponent. Following the rules formulated above, we first make a guess about the form of the polynomial part of a particular solution $y_p(x)$ and then multiply it by the exponent:$$y_p(x) = (Ax+B)e^{2x}$$To obtain the coefficients, we substitute this "half-baked" solution into the original equation:$$y_p' = Ae^{2x} + 2e^{2x}(Ax+B) \\ y_p'' = 2Ae^{2x} + 2(Ae^{2x} + 2e^{2x}(Ax+B)) = 4Ae^{2x} + 4e^{2x}(Ax+B)\\ y_p'' -2 y_p +10y_p = xe^{2x}\\ 4Ae^{2x} + 4e^{2x}(Ax+B) - 2Ae^{2x} - 4e^{2x}(Ax+B) + 10e^{2x}(Ax+B) = xe^{2x}\\ 4A + 4Ax+4B-2A-4Ax-4B+10Ax+10B = x\\ 10Ax + 2A+10B = 1\cdot x+0\\$$To obtain the coefficients $A$ and $B$, we need to equate the coefficients of corresponding powers of x. This means,$$\begin{cases} 10A = 1\\ 2A+10B = 0 \end{cases} \implies \begin{cases} A = \frac{1}{10}\\ B = -\frac{1}{50} \end{cases}$$So finally, our particular solution is$$y_p(x) = (\frac{1}{10}x - \frac{1}{50})e^{2x}$$
3. Step 3: gathering both parts into the general solution:$$y(x) = y_0(x)+y_p(x) = e^{-x}(C_1cos3x + C_2sin3x) + (\frac{1}{10}x - \frac{1}{50})e^{2x}, \\ C_1, C_2 = const$$

Congratulations! We've solved a second-order nonhomogeneous ODE.

Conclusion

Differential equations of higher orders are more complicated than linear ones, but the same algorithms may be applied to solve them. First, we need to find the solution to the homogeneous equation by solving the accompanying characteristic equation. Its roots completely determine the form of the solution. Then you need to find a particular solution to the nonhomogeneous equation – there are several heuristics for this. Finally, the two parts must be combined to obtain a general solution.