Skip to content
Snippets Groups Projects

Higher order linear differential equations

In the previous lecture, we focused on first order linear differential equations as well as systems of such equations. In this lecture we switch focus to DE's which involve higher derivatives of the function we would like to solve for. To f`%acilitate this change we are going to change notation. In the previous lecture we wrote differential equations for

x(t)x(t)
. In this lecture we will write DE's of
y(x)y(x)
, where
yy
is an unknown function and
xx
is the independent variable. For this purpose we make the following definitions,

y=dydx, y=d2ydx2, , y(n)=dnydxn.y' = \frac{dy}{dx}, \ y'' = \frac{d^2 y}{dx^2}, \ \cdots, \ y^{(n)} = \frac{d^n y}{dx^n}.

In the new notation, a linear

nn
-th order differential equation with constant coefficients reads

y(n)+an1y(n1)++a1y+a0y=0.y^{(n)} + a_{n-1} y^{(n-1)} + \cdots + a_1 y' + a_0 y = 0.

!!! info "Linear combination of solutions are still solutions"

Note that as was the case for first order linear DE's, the propery of 
linearity once again means that if $y_{1}(x)$ and $y_{2}(x)$ are both 
solutions, and $a$ and $b$ are constants, 

$$a y_{1}(x) + b y_{2}(x)$$

then linear combination of the solutions is also a solution.

In order to solve a higher order linear DE we will present a trick that makes it possible to map the problem of solving a single

nn
-th order linear DE into a related problem of solving a system of
nn
first order linear DE's.

To begin, define:

y1=y, y2=y, , yn=y(n1).y_{1} = y, \ y_{2} = y', \ \cdots, \ y_{n} = y^{(n-1)}.

Then, the differential equation can be re-written as

y1=y2y_1 ' = y_2
y2=y3y_2 ' = y_3
\vdots
yn1=yny_{n-1} ' = y_{n}
yn=a0y1a1y2an1yn.y_{n} ' = - a_{0} y_{1} - a_{1} y_{2} - \cdots - a_{n-1} y_{n}.

Notice that together these

nn
equations form a linear first order system, the first
n1n-1
equations of which are trivial. Note that this trick can be used to reduce any system of
nn
-th order linear DE's to a larger system of first order linear DE's.

Since we have discussed already the method of solution for first order linear systems, we will outline the general solution to this system. As before, the general solution will be the linear combination of

nn
linearly independent solutions
fi(x)f_{i}(x)
,
iϵ{1,,n}i \epsilon \{1, \cdots, n \}
, which make up a basis for the solution space. That is the general solution has the form

y(x)=c1f1(x)+c2f2(x)++cnfn(x).y(x) = c_1 f_1 (x) + c_2 f_2 (x) + \cdots + c_n f_{n}(x).

To check that the

nn
solutions form a basis, it is sufficient to verify

det[f1(x)fn(x)f1(x)fn(x)f1(n1)(x)fn(n1)(x)]0. det \begin{bmatrix} f_1(x) & \cdots & f_{n}(x) \\ f_1 ' (x) & \cdots & f_{n}'(x) \\ \vdots & \vdots & \vdots \\ f^{(n-1)}_{1} (x) & \cdots & f^{(n-1)}_{n} (x) \\ \end{bmatrix} \neq 0.

The determinant in the preceding line is called the Wronski determinant. In particular, to determine solutions, we need to find the eigenvalues of

A=[01000010a0a1a2an1].A = \begin{bmatrix} 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ -a_0 & -a_1 & -a_2 & \cdots & -a_{n-1} \\ \end{bmatrix}.

It is possible to show that

det(AλI)=P(λ),det(A - \lambda I) = -P(\lambda),

in which

P(λ)P(\lambda)
is the characteristic polynomial of the system matrix
AA
,

P(λ)=λn+an1λn1++a0.P(\lambda) = \lambda^n + a_{n-1} \lambda^{n-1} + \cdots + a_0.

As we demonstrate below, the proof relies on the co-factor expansion technique for calculating a determinant.

det(AλI)=[λ1000λ10a0a1a2an1+λ]- det(A - \lambda I) = \begin{bmatrix} \lambda & -1 & 0 & \cdots & 0 \\ 0 & \lambda & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ a_0 & a_1 & a_2 & \cdots & a_{n-1} + \lambda \\ \end{bmatrix}
det(AλI)=λdet[λ1000λ10a1a2a3an1+λ]+(1)n+1a0det[1000λ10cdots000λ1]- det(A - \lambda I) = \lambda det \begin{bmatrix} \lambda & -1 & 0 & \cdots & 0 \\ 0 & \lambda & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ a_1 & a_2 & a_3 & \cdots & a_{n-1} + \lambda \\ \end{bmatrix} + (-1)^{n+1}a_0 det \begin{bmatrix} -1 & 0 & 0 & \cdots & 0 \\ \lambda & -1 & 0 & cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 \\ \end{bmatrix}
det(AλI)=λdet[λ1000λ10a1a2a3an1+λ]+(1)n+1a0(1)n1- det(A - \lambda I) = \lambda det \begin{bmatrix} \lambda & -1 & 0 & \cdots & 0 \\ 0 & \lambda & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ a_1 & a_2 & a_3 & \cdots & a_{n-1} + \lambda \\ \end{bmatrix} + (-1)^{n+1} a_0 (-1)^{n-1}
det(AλI)=λdet[λ1000λ10a1a2a3an1+λ]+a0- det(A - \lambda I) = \lambda det \begin{bmatrix} \lambda & -1 & 0 & \cdots & 0 \\ 0 & \lambda & -1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ a_1 & a_2 & a_3 & \cdots & a_{n-1} + \lambda \\ \end{bmatrix} + a_0
det(AλI)=λ(λ(λ+a2)+a1)+a0- det(A - \lambda I) = \lambda (\lambda (\lambda \cdots + a_2) + a_1) + a_0
det(AλI)=P(λ).- det(A - \lambda I) = P(\lambda).

In the second last line of the proof we indicated that the method of co-factor expansion demonstrated is repeated an additional

n2n-2
times. This completes the proof.

With the characteristic polynomial, it is possible to write the differential equation as

P(ddx)y(x)=0.P(\frac{d}{dx})y(x) = 0.

To determine solutions, we need to find

λi\lambda_i
such that
P(λi)=0P(\lambda_i) = 0
. By the fundamental theorem of algebra, we know that
P(λ)P(\lambda)
can be written as

P(λ)=Σk=1l(λλk)mk.P(\lambda) = \overset{l}{\underset{k=1}{\Sigma}} (\lambda - \lambda_k)^{m_k}.

In the previous equation

λk\lambda_k
are the k roots of the equations, and
mkm_k
is the multiplicity of each root. Note that the multiplicities satisfy
Σk=1lmk=n\overset{l}{\underset{k=1}{\Sigma}} m_k = n
.

If the multiplicity of each eigenvalue is one, then solutions which form the basis are then given as:

f(x)=eλ1x, eλ2x, , eλnx.f(x) = e^{\lambda_1 x}, \ e^{\lambda_2 x}, \ \cdots, \ e^{\lambda_n x}.

If there are eigenvalues with multiplicity greater than one, the the solutions which form the basis are given as

f(x)=eλ1x, xeλ1x, , xm11eλ1x, etc.f(x) = e^{\lambda_1 x}, \ x e^{\lambda_1 x} , \ \cdots, \ x^{m_{1}-1} e^{\lambda_1 x}, \ etc.

---ADD PROOF HERE---

!!! check "Example: Second order homogeneous linear DE with constant coefficients"

Consider the equation 

$$y'' + Ey = 0.$$ 

The characteristic polynomial of this equation is 

$$P(\lambda) = \lambda^2 + E.$$

There are three cases for the possible solutions, depending upon the value 
of E.

**Case 1: $E>0$**
For ease of notation, define $E=k^2$ for some constant $k$. The 
characteristic polynomial can then be factored as

$$P(\lambda) = (\lambda+ i k)(\lambda - i k). $$

Following our formulation for the solution, the two basis functions for the 
solution space are 

$$f_1(x) = e^{i k x}, \ f_2=e^{- i k x}.$$

Alternatively, the trigonometric functions can serve as basis functions, 
since they are linear combinations of $f_1$ and $f_2$ which remain linearly
independent,

$$\tilde{f_1}(x)=cos(kx), \tilde{f_2}(x)=sin{kx}.$$

**Case 2: $E<0$**
This time, define $E=-k^2$, for constant $k$. The characteristic polynomial 
can then be factored as 

$$P(\lambda) = (\lambda+ k)(\lambda -  k).$$

The two basis functions for this solution are then 

$$f_1(x)=e^{k x}, \ f_2(x) = e^{-k x}.$$

**Case 3: $E=0$**
In this case, there is a repeated eigenvalue (equal to $0$), since the 
characteristic polynomial reads

$$P(\lambda) = (\lambda-0)^2.$$

Hence the basis functions for the solution space read 

$$f_1(x)=e^{0 x} = 1, \ f_{2}(x) = x e^{0 x} = x. $$

Partial differential equations

A partial differential equation (PDE) is an equation involving a function of two or more indepenedent variables and derivatives of said function. These equations are classified similarly to ordinary differential equations (the subject of our earlier study). For example, they are called linear if no terms such as

y(x,t)xdy(x,t)t or\frac{\partial y(x,t)}{\partial x} \cdot \frac{d y(x,t)}{\partial t} \ or
2y(x,t)x2y(x,t)\frac{\partial^2 y(x,t)}{\partial x^2} y(x,t)

occur. A PDE can be classified as

nn
-th order accorind to the highest derivative order of either variable occuring in the equation. For example, the equation

3f(x,y)x3+f(x,t)t=5\frac{\partial^3 f(x,y)}{\partial x^3} + \frac{\partial f(x,t)}{\partial t} = 5

is a 3-rd order equation because of the third derivative with respect to x in the equation.

To begin, we demonstrate that PDE's are of fundamental importance in physics, especially in quantum physics. In particular, the Schrödinger equation, which is of central importance in quantum physics is a partial differential equation with respect to time and space. This equation is very important because it describes the evolution in time and space of the entire description of a quantum system \psi(x,t), which is known as the wavefunction.

For a free particle in one dimension, the Schrödinger equation is

i \hbar \frac{\partial \psi(x,t)}{\partial t} = - \frac{\hbar^2}{2m} \frac{\partial^2 \psi(x,t)}{\partial x^2}.

When we studied ODEs, an initial condition was necessary in order to fully specify a solution. Similarly, in the study of PDEs an initial condition is required but now boundary conditions are also required. Going back to the intuitive discussion from the lecture on ODEs, each of these conditions is necassary in order to specify an integration constant that occurs in solving the equation. In partial differential equations at least one such constant will arise from the time derivative and likewise at least one from the spatial derivative.

For the Schrödinger equation, we could supply the initial conditions

\psi(x,0)= \psi_{0}(x) \ & \ \psi(0,t) = \psi{t, L} = 0.

This particular set of boundary conditions corresponds to a particle in a box, a situation which is used as the base model for many derivations in quantum physics.

Another example of a partial differential equation common in physics is the Laplace equation

\frac{\partial^2 \phi(x,y)}{\partial x^2}+\frac{\partial^2 \phi(x,y)}{\partial y^2}=0.

In quantum physics Laplace's equation is important for the study of the hydrogen atom. In three dimensions and using spherical coordinates, the solutions to Laplace's equation are special functions called spherical harmonics. In the context of the hydrogen atom, these functions describe the wave function of the system and a unique spherical harmonic function corresponds to each distinct set of quantum numbers.

In the study of PDEs there is not a comprehensive overall treatment to the same extent as there is for ODEs. There are several techniques which can be applied to solving these equations, but the choice of technique must be tailored to the equation at hand. Hence we focus on some specific examples that are common in physics.

Separation of variables

Let us focus on the one dimensional Schrödinger equation of a free particle

i \hbar \frac{\partial \psi(x,t)}{\partial t} = - \frac{\hbar^2}{2m} \frac{\partial^2 \psi(x,t)}{\partial x^2}.

To attempt a solution, we will make a separation ansatz,

\psi(x,t)=\phi(x) f(t).

!!! info "Separation ansatz" The separation ansatz is a restrictive ansatz, not a fully general one. In general, for such a treatment to be valid an equation and the boundary conditions given with it have to fulfill certain properties. In this course however you will only be asked to use this technique when it is suitable.

Substituting the separation ansatz into the PDE,

i \hbar \frac{\partial \phi(x)f(t)}{\partial t} = - \frac{\hbar^2}{2m} \frac{\partial^2 \phi(x)f(t)}{\partial x^2} i \hbar \dot{f}(t) \phi(x) = - \frac{\hbar^2}{2m} \phi''(x)f(t).

Notice that in the above equation the derivatives on f and \phi can each be written as ordinary derivatives, \dot{f}=\frac{df(t)}{dt}, \phi''(x)=\frac{d^2 \phi}{dx^2}. This is so because each is only a function of one variable.

Next, divide both sides of the equation through by \psi(x,t)=\phi(x) f(t),

i \hbar \frac{\dot{f}(t)}{f(t)} = - \frac{\hbar^2}{2m} \frac{\phi''(x)}{\phi(x)} = constant := \lambda.

In the previous line we concluded that each part of the equation must be equal to a constant, which we defined as \lambda. This follows because the left hand side of the equation only has a dependence on the spatial coordinate x, whereas the right hand side only has dependence on the time coordinate t. If we have two functions a(x) and b(t) such that a(x)=b(t) \ \forall x, \ t \ \epsilon \mathbb{R}, then a(x)=b(t)=const.

The constant we defined, lambda, is called a separation constant. With it, we can break the spatial and time dependent parts of the equation into two separate equations,

i \hbar \dot{f}(t) = \lambda f(t)

-\frac{\hbar^2}{2m} \phi''(x) = \lambda \phi(x) .

To summarize, this process has broken one partial differential equation into two ordinary differential equations of different variables. In order to do this, we needed to introduce a separation constant, which remains to be determined.

Boundary and eigenvalue problems

Continuing on with the Schrödinger equation example from the previous section, let us focus on

-\frac{\hbar^2}{2m} \phi''(x) = \lambda \phi(x), \phi(0)=\phi(L)=0.

This has the form of an eigenvalue equation, in which \lambda is the eigenvalue, - \frac{\hbar^2}{2m} \frac{d^2}{dx^2}[\cdot] is the linear operator and \phi(x) is the eigenfunction.

Notice that when stating the ordinary differential equation, it is specified along with it's boundary conditions. Note that in contrast to an initial value problem, a boundary value problem does not always have a solution. For example, in the figure below, regardless of the initial slope, the curves never reach 0 when x=L.

For boundary value problems like this, there are only solutions for particular eigenvalues \lambda. Coming back to the example, it turns out that solutions only exist for \lambda>0 --this can be shown quickly, feel free to try it! Define for simplicity k^2:= \frac{2m \lambda}{\hbar^2}. The equation then reads

\phi''(x)+k^2 \phi(x)=0.

Two linearly independent solutions to this equation are

\phi_{1}(x)=sin(k x), \ \phi_{2}(x) = cos(k x).

The solution to this homogeneous equation is then

\phi(x)=c_1 \phi_1(x)+c_2 \phi_2(x).

The eigenvalue, \lambda as well as one of the constant coefficients can be determined using the boundary conditions.

\phi(0)=0 \ \Rightarrow \ \phi(x)=c_1 sin(k x), \ c_2=0.

\phi(L)=0 \ \Rightarrow \ 0=c_1 sin(k L) .

In turn, using the properties of the sin(\cdot) function, it is now possible to find the allowed values of k and hence also \lambda. The previous equation implies,

k L = n \pi, \ n \ \epsilon \ \mathbb{N}

\lambda_n = \big{(}\frac{n \pi \hbar}{L} \big{)}^2.

The values \lambda_n are the eigenvalues. Now that we have determined \lambda, it enters into the time equation, i \hbar \dot{f}(t) = \lambda f(t) only as a constant. We can hence simply solve,

\dot{f}(t) = -i \frac{\lambda}{\hbar} f(t)

f(t) = A e^{\frac{-i \lambda t}{\hbar}}.

In the previous equation, the coefficient A can be determined if the original PDE was supplied with an initial condition.

Putting the solutions to the two ODEs together and redefining \tilde{A}=A \cdot c_1, we arrive at the solutions for theb PDE,

\psi_n(x,t) = \tilde{A}_n e^{-i \frac{\lambda_n t}{\hbar}} sin(\frac{n \pi x}{L}).

Notice that there is one solution \psi_{n}(x,t) for each natural number n. These are still very special solutions. We will begin discussing next how to obtain the general solution in our example.

Self-adjoint differential equations: Connection to Hilbert spaces!

As we hinted was possible earlier, let us re-write the previous equation by defining a linear operator, L, acting on the space of functions which satisfy \phi(0)=\phi(L)=0:

L[\cdot]:= \frac{- \hbar^2}{2m} \frac{d^2}{dx^2}[\cdot].

Then, the ODE can be writted as

L[\phi]=\lambda \phi.

This equation looks exactly like, and turns out to be, an eigenvalue equation!

!!! info "Connecting function spaces to Hilbert spaces"

Recall that a space of functions can be transformed into a Hilbert space by 
equipping it with a inner product,

$$\langle f, g \rangle = \int^{L}_{0} dx f*(x) g(x) $$

Use of this inner product also has utility in demonstrating that particular 
operators are *Hermitian*. The term hermitian is precisely defined below.
Of considerable interest is that hermition operators have a set of nice 
properties including all real eigenvalues and orthonormal eigenfunctions. 

The nicest type of operators for many practical purposes are hermitian operators. In quantum physics for example, any physical operator must be hermitian. Denote a hilbert space \mathcal{H}. An opertor H: \mathcal{H} \mapsto \mathcal{H} is said to be hermitian if it satisfies

\langle f, H g \rangle = \langle H f, g \rangle \ \forall \ f, \ g \ \epsilon \ \mathcal{H}.

Now, we would like to investigate whether the operator we have been working with, L satisfies the criterion of being hermitian over the function space \phi(0)=\phi(L)=0 equipped with the above defined inner product (i.e. it is a Hilbert space). Denote this Hilbert space \mathcal{H}_{0} and consider let f, \ g \ \epsilon \ \mathcal{H}_0 denote two functions from the Hilbert space. Then, we can investigate

\langle f, L g \rangle = \frac{- \hbar^2}{2m} \int^{L}_{0} dx f*(x) \frac{d^2}{dx^2}g(x).

As a first step, it is possible to do integration by parts in the integral,

\langle f, L g \rangle = \frac{+ \hbar^2}{2m} ( \int^{L}_{0} dx \frac{d f*}{dx} \frac{d g}{dx} - [f*(x)\frac{d g}{dx}] \big{|}^{L}_{0} )

The boundary term vansishes, due to the boundary conditions f(0)=f(L)=0, which directly imply f*(0)=f*(L)=0. Now, intergrate by parts a second time

\langle f, L g \rangle = \frac{- \hbar^2}{2m} (\int^{L}_{0} dx \frac{d^2 f*}{dx^2} g(x) - [\frac{d f*}{dx} g(x)] \big{|}^{L}_{0} ).

As before, the boundary term vanishes, due to the boundary conditions g(0)=g(L)=0. Upon cancelling the boundary term however, the expression on the right hand side, contained in the integral is simply \langle L f, g \rangle. Therefore,

\langle f, L g \rangle=\langle L f, g \rangle.

We have demonstrated that L is a hermitian operator on the space \mathcal{H}_0. As a hermitian operator, L has the property that it's eigenfunctions form an orthonormal basis for the space \mathcal{H}_0. Hence it is possible to expand any function f \ \epsilon \ \mathcal{H}_0 in terms of the eigenfunctions of L.

!!! info "Connection to quantum states"

Recall that q quantum state $\ket{\phi}$ can be written in an orthonormal 
basis $\{ \ket{u_n} \}$ as 
$$\ket{\phi} = \underset{n}{\Sigma} \bra{u_n} \ket{\phi} \ket{u_n}.$$ 

In terms of hermitian operators and their eigenfunctions, the eigenfunctions
play the role of the orthonormal basis. In reference to our running example,
the 1D Schrödinger equation of a free particle, the eigenfunctions 
$sin(\frac{n \pi x}{L})$ play the role of the basis functions $\ket{u_n}$.

To close our running example, consider the initial condition \psi(x,o) = \psi_{0}(x). Since the eigenfunctions sin(\frac{n \pi x}{L}) form a basis, we can now write the general solution to the problem as

\psi(x,t) = \overset{\infinity}{\underset{n}{\Sigma}} c_n e^{-i \frac{\lambda_n t}{\hbar}} sin(\frac{n \pi x}{L}),

where in the above we have defined the coefficients using the Fourier coefficient,

c_n:= \int^{L}_{0} dx sin(\frac{n \pi x}{L}) \psi_{0}(x).

General recipie for seperable PDEs

  1. Make the separation ansatz to obtain separate ordinary differential equations.
  2. Choose which euation to treat as the eigenvalue equation. This will depend upon the boundary conditions. Additionally, verify that the linear differential operator L in the eigenvalue equation is hermitian. 3.Solve the eigenvalue equation. Substitute the eigenvalues into the other equations and solve those too.
  3. Use the orthonormal basis functions to write down the solution corresponding to the specified initial and boundary conditions.

One natural question is what if the operator L from setp 2 is not hermitian? It is possible to try and make it hermitian by working on a Hilbert space equipped with a different inner product. This means one can consider modifications to the definition of \langle \cdot, \cdot \rangle such that L is hermitian with respect to the modified inner product. This type of technique falls under the umbrella of Sturm-Liouville Theory, which forms the foundation for a lot of the analysis that can be done analytically on PDEs.

Another question is of course what if the equation is not separable? One possible approach is to try working in a different coordinate system. There are a few more analytic techniques available, however in many situations it becomes necessary to work with numerical methods of solution.

Problems

  1. [:grinning:] Which of the following equations for y(x) is linear?

    (a) y''' - y'' + x cos(x) y' + y - 1 = 0
    
    (b) y''' + 4 x y' - cos(x) y = 0
    
    (c) y'' + y y' = 0
    
    (d) y'' + e^x y' - x y = 0
  2. [:grinning:] Find the general solution to the equation

    $$y'' - 4 y' + 4 y = 0. $$
    
    Show explicitly by computing the Wronski determinant that the 
    basis for the solution space is actually linearly independent. 
  3. [:grinning:] Find the general solution to the equation

    $$y''' - y'' + y' - y = 0.$$
    
    Then find the solution to the initial conditions $y''(0) =0$, $y'(0)=1$, $y(0)=0$. 
  4. [:smirk:] Take the Laplace equation in 2D:

    $$\frac{\partial^2 \phi(x,y)}{\partial x^2} + \frac{\partial^2 \phi(x,y)}{\partial y^2} = 0.$$
    
    (a) Make a separation ansatz $\phi(x,y) = f(x)g(y)$ and write 
      down the resulting ordinary differential equations.
    
    (b) Now assume that the boundary conditions $\phi(0,y) = \phi(L,y) =0$ for 
      all y, i.e.  f(0)=f(L)=0. Find all solutions $f(x)$ and the corresponding 
      eigenvalues.
    
    (c) Finally, for each eigenvalue, find the general solution $g(y)$ for this
      eigenvalue. Combine this with all solutions $f(x)$ to write down the general
      solution (we know from the lecture that the operator $\frac{d^2}{dx^2}$ is 
      hermitian - you can thus directly assume that the solutions form an orthogonal
      basis). 
  5. [:smirk:] Take the partial differential equation

    $$\frac{\partial h(x,y)}{\partial x} + x \frac{\partial h(x,y)}{\partial y} = 0. $$
    
    Try to make a separation ansatz $h(x,y)=f(x)g(y)$. What do you observe?
  6. [:sweat:] Bonus question - this kind of question will not be asked in the exam

    We consider the Hilbert space of functions $f(x)$ defined for $x \ \epsilon \ [0,L]$
    with $f(0)=f(L)=0$. 
    
    Which of the following operators on this space is hermitian?
    
    (a) Lf = A(x) \frac{d^2 f}{dx^2}
    
    (b) Lf = \frac{d}{dx} \big{()} A(x) \frac{df}{dx} \big{)}