Sturm–Liouville theory - show that the linear operator is self-adjoint.

The goal of this post is to document my derivation that the linear operator for the Sturm Liouville theory is actually self-adjoint

The differential equation is

$ \frac{d}{dx}(p(x) \frac{d}{dx}y) + (\lambda w(x) + q(x))y = 0 $

Rearranging, we can write this as an eigenvalue problem as follow:

$ \frac{d}{dx}(p \frac{d}{dx}y) + qy = -\lambda w y $

$ \frac{1}{w}(\frac{d}{dx}(p \frac{d}{dx}y) + qy) = -\lambda y $

Define the operator $ L = \frac{1}{w}(\frac{d}{dx}(p \frac{d}{dx}) + q) $, it is trivial to verify the operator is linear.

Consider the function space with inner product defined as $ <a, b> = \int a b w dx $, we wanted to verify the operator $ L $ is self-adjoint in this space, in other words, $ <As, t> = <s, At> $.

$ \begin{eqnarray*} <As, t> &=& \int \frac{1}{w}(\frac{d}{dx}(p \frac{d}{dx}s) + qs) t w dx \\ &=& \int (\frac{d}{dx}(p \frac{d}{dx}s) + qs) t dx \\ &=& \int \frac{d}{dx}(p \frac{d}{dx}s)t dx + \int qst dx \end{eqnarray*} $.

Similarly we can write the other side as

$ \begin{eqnarray*} <s, At> &=& \int \frac{1}{w}(\frac{d}{dx}(p \frac{d}{dx}t) + qt) s w dx \\ &=& \int (\frac{d}{dx}(p \frac{d}{dx}t) + qt) s dx \\ &=& \int \frac{d}{dx}(p \frac{d}{dx}t)s dx + \int qts dx \end{eqnarray*} $.

So all we needed to prove is $ \int (ps')'t dx = \int (pt')'s dx $, we do that using integration by part:

$ \begin{eqnarray*} & & \int (ps')'t dx \\ &=& \int t d(ps') dx \\ &=& ps't - \int ps' dt \\ &=& ps't - \int ps't' dx \\ &=& ps't - \int pt' ds \\ &=& ps't -pt's + \int s d(pt') \\ &=& ps't -pt's + \int (pt')' s dx \\ \end{eqnarray*} $.

So we need to show $ ps't - pt's = 0 $. This should come from the boundary conditions, which is not specified here. The boundary of a Sturm Liouville problem should satisfy this. For example, $ p(x) = 0 $ for the Legendre's problem at the boundary.

The use of a known solution to find another

Problem:

$ y'' - \frac{1}{2x}y' - \frac{1}{x^2}y = 0 $

Solution:

A good guess would be $ y = x^{k+2} $ with an unknown $ k $, the $ k + 2 $ is really just to make the algebra simpler.

$ \begin{eqnarray*} y'' - \frac{1}{2x}y' - \frac{1}{x^2}y &=& 0 \\ (k+2)(k+1)x^k - \frac{1}{2x}(k+2)(x^{k+1}) - \frac{1}{x^2}x^{k+2} &=& 0 \\ (k+2)(k+1) - \frac{k+2}{2} - 1 &=& 0 \\ 2(k+2)(k+1) - (k+2) - 2 &=& 0 \\ 2k^2 + 5k &=& 0 \\ k &=& 0 \text{ or } \frac{-5}{2} \\ \end{eqnarray*} $

So the answer is $ y = Ax^2 + \frac{B}{\sqrt{x}} $

For the sake of exercise, let's assume we do not know the second independent answer $ y = \frac{1}{\sqrt{x}} $. Can we use the answer $ y = x^2 $ to derive it?

Here is a technique for using an existing solution for another one. Assume an answer for $ y'' + p(x)y' + q(x)y = 0 $ is $ y_1(x) $, we let $ y_2(x) = v(x)y_1(x) $ and substitute this back to the equation, we get:

$ \begin{eqnarray*} y'' - py' - qy &=& 0 \\ (v''y_1 + 2v'y_1' + vy_1'') + p(v'y_1 + vy_1') + q(vy_1) &=& 0 \\ v''y_1 + 2v'y_1' + p(v'y_1) + vy_1'' + p(vy_1') + q(vy_1) &=& 0 \\ v''y_1 + 2v'y_1' + p(v'y_1) &=& 0 \\ \frac{v''}{v'} + 2\frac{y_1'}{y_1} + p &=& 0 \\ \frac{v''}{v'} &=& - 2\frac{y_1'}{y_1} - p \\ \log v' &=& -2 \log y_1 - \int p dx \\ v' &=& e^{-2 \log y_1 - \int p dx} \\ &=& \frac{1}{y_1^2}e^{-\int p dx} \\ v &=& \int \frac{1}{y_1^2}e^{-\int p dx} dx \end{eqnarray*} $

So we will need to do is to compute the two integral for our special problem.

$ \begin{eqnarray*} v &=& \int \frac{1}{y_1^2}e^{-\int p dx} dx \\ &=& \int \frac{1}{x^4}e^{-\int \frac{-1}{2x} dx} dx\\ &=& \int \frac{1}{x^4}e^{-(\frac{-1}{2}\log x)} dx \\ &=& \int \frac{1}{x^4}e^{\frac{1}{2}\log x} dx \\ &=& \int \frac{1}{x^4}e^{\log \sqrt{x}} dx \\ &=& \int \frac{1}{x^4}\sqrt{x} dx \\ &=& \int x^{-3.5} dx \\ &=& \frac{x^{-2.5}}{2.5} \\ y_2 &=& vy_1 \\ &=& \frac{x^{-2.5}}{2.5}x^2 \\ &=& \frac{x^{-0.5}}{2.5} \end{eqnarray*} $

As constant can be ignored, now we come back to full circle. In general, I guess the integrals are going to be hard to evaluate. This problem seems rather tailored for the purpose of practicing.

The method of variation of parameter (II)

Problem:

$ y'' - 3y' + 2y = \cos x $

Solution:

The solution to the homogeneous equation $ w'' - 3w' + 2w = 0 $ is obviously $ w_1 = e^x $ and $ w_2 = e^{2x} $.

Next, we solve the equation pair:

$ v_1' e^x + v_2' e^{2x} = 0 $.
$ v_1' e^x + 2 v_2' e^{2x} = \cos x $.

The second equation subtracts the first equation gives:

$ v_2' e^{2x} = \cos x $

Therefore $ v_2' = \cos x e^{-2x} $.

This time I do the integral myself using Euler's relation!

Consider the integral $ \int{e^{ix}e^{-2x}dx} = \int{(\cos x + i\sin x)e^{-2x}dx} $. Therefore we can equate the real part to get the result, and the first integral is particularly easy.

$ \int{e^{ix}e^{-2x}dx} = \frac{1}{-2 + i}e^{(-2 + i)x} $.

Next we find the real part of it by rationalizing and Euler's relation.

$ \begin{eqnarray*} & & \frac{1}{-2 + i}e^{(-2 + i)x} \\ &=& \frac{(-2 - i)}{(-2 + i)(-2 - i)}e^{(-2 + i)x} \\ &=& \frac{(-2 - i)}{5}e^{(-2 + i)x} \\ &=& \frac{(-2 - i)}{5}e^{-2x}e^{ix} \\ &=& \frac{(-2 - i)}{5}e^{-2x}(\cos x + i\sin x) \\ \end{eqnarray*} $

Therefore, we have

$ \begin{eqnarray*} & & v_2(x) \\ &=& \int \cos x e^{-2x} dx \\ &=& e^{-2x}(\frac{-2}{5}\cos x + \frac{1}{5}\sin x) \end{eqnarray*} $

Similarly, we solve $ v_1'(x) $ by subtracting the second equation by twice the first one, that gives $ v_1'(x) = -\cos(x)e^{-2x} $. The integral is very similar and we get $ \begin{eqnarray*} & & v_1(x) \\ &=& \int -\cos x e^{-x} dx \\ &=& e^{-x}(\frac{1}{2}\cos x - \frac{1}{2}\sin x) \end{eqnarray*} $

Therefore the answer is $ (\frac{-2}{5} + \frac{1}{2})\cos x + (\frac{1}{5} - \frac{1}{2})\sin x + Ae^x + Be^{2x} = \frac{1}{10}\cos x - \frac{3}{10}\sin x + Ae^x + Be^{2x} $.

The method of variation of parameter (I)

Problem:

$ y'' - 2y = x + 1 $.

Solution:

Let's discuss the method in general first, consider the differential equation:

$ y'' + p(x)y' + q(x)y = r(x) $.

Assume we can solve the homogeneous equation

$ w'' + p(x)w' + q(x)w = 0 $

and get the solution the two independent solutions $ w(x) = w_1(x) $ and $ w(x) = w_2(x) $.

Now we assume the solution to the original equation takes the form:

$ y(x) = v_1(x)w_1(x) + v_2(x)w_2(x) $.

Computing the first derivative, we get:

$ y'(x) = v_1'(x)w_1(x) + v_1(x)w_1'(x) + v_2'(x)w_2(x) + v_2(x)w_2'(x) $.

As a wishful thinking, let's assume $ v_1'(x)w_1(x) + v_2'(x)w_2(x) = 0 $. This is one of the equation that govern $ v_1'(x) $ and $ v_2'(x) $.

Next, we compute the second derivative, we get:

$ y'(x) = v_1(x)w_1'(x) + v_2(x)w_2'(x) $.
$ y''(x) = v_1'(x)w_1'(x) + v_1(x)w_1''(x) + v_2'(x)w_2'(x)  + v_2(x)w_2''(x) $.

Putting these derivative back to the differential equation, we get:

$ (v_1'(x)w_1'(x) + v_1(x)w_1''(x) + v_2'(x)w_2'(x)  + v_2(x)w_2''(x)) + p(x)(v_1(x)w_1'(x) + v_2(x)w_2'(x)) + q(x)(v_1(x)w_1(x) + v_2(x)w_2(x)) = r(x) $.

A very long equation, but if we carefully group the terms we see this:

$ v_1'(x)w_1'(x) + v_2'(x)w_2'(x) + (v_1(x)w_1''(x) + p(x)v_1(x)w_1'(x) + q(x)v_1(x)w_1(x)) + (v_2(x)w_2''(x) + p(x)v_2(x)w_2'(x) + q(x)v_2(x)w_2(x)) = r(x) $.

Note that the two bracketed term is actually 0 because $ w_1(x) $ and $ w_2(x) $ is a solution to the homogeneous equation. So it simplifies to

$ v_1'(x)w_1'(x) + v_2'(x)w_2'(x) = r(x) $.

That will be our second equation that govern $ v_1'(x) $ and $ v_2'(x) $, with this two equations, we can solve the $ v_1(x) $ and $ v_2(x) $ and therefore the overall differential equation.

So much for the theory, now we solve the equation $ y'' - 2y = x + 1 $.

The homogeneous equation $ w'' - 2w = 0 $ is pretty easy to solve, it is simply $ w_1(x) = e^{\sqrt{2}x} $ and $ w_2(x) = e^{-\sqrt{2}x} $.

Next, we use the two equations to solve for the $ v $ functions

$ v_1'(x)e^{\sqrt{2}x} + v_2'(x)e^{-\sqrt{2}x} = 0 $
$ \sqrt{2}v_1'(x)e^{\sqrt{2}x} + (-\sqrt{2})v_2'(x)e^{-\sqrt{2}x} = x+1 $

Cancelling the $ \sqrt{2} $ factor from the second equation we get

$ v_1'(x)e^{\sqrt{2}x} - v_2'(x)e^{-\sqrt{2}x} = \frac{x + 1}{\sqrt{2}} $

Now the functions are trivial to solve $ v_1'(x) = \frac{x+1}{2\sqrt{2}}e^{-\sqrt{2}x} $, $ v_2'(x) = \frac{-x-1}{2\sqrt{2}}e^{\sqrt{2}x} $.

Next we integrate to find $ v_1(x) $ and $ v_2(x) $.

$ v_1(x) = \frac{-1}{8}(2x + \sqrt{2} + 2) e^{-\sqrt{2}x} $
$ v_2(x) = \frac{-1}{8}(2x - \sqrt{2} + 2) e^{\sqrt{2}x} $

Putting these back to the differential equation, we finally get

$ y(x) = v_1(x)w_1(x) + v_2(x)w_2(x) = \frac{-1}{8}(2x + \sqrt{2} + 2) +  \frac{-1}{8}(2x - \sqrt{2} + 2) = \frac{-1}{2}(x + 1) $.

And the final answer is, of course, $ y(x) = \frac{-1}{2}(x + 1) + Ae^{\sqrt{2}x} + Be^{-\sqrt{2}x} $.

To be honest to myself, I cheated. I didn't do all these integration myself. I used these simple MATLAB code. It is really troublesome to use this method.

clc
w1 = exp(sqrt(2) * x);
w2 = exp(-sqrt(2) * x);
v1d = (x + 1)/sqrt(2) /2 / w1;
v2d = -(x + 1)/sqrt(2) /2 / w2;
g = simplify(int(v1d) * w1 + int(v2d) * w2)
diff(diff(g)) - 2 * g

The method of undetermined coefficients (II)

Problem:

y'' + y' = x

Solution:

Again, we solve the homogeneous equation $ z'' + z' = 0 $, the solution is $ z = Ae^{-x} + B $.

Next we guess the particular solution to $ w'' + w' = x $ is a quadratic polynomial $ w(x) = Px^2 + Qx + R $, the derivatives are

$ w'(x) = 2Px + Q $
$ w''(x) = 2P $.

Therefore we have $ 2P + 2Px + Q = x $, so we get $ P = \frac{1}{2} $ and $ Q = -1 $.

The final answer is therefore $ y = Ae^{-x} + \frac{1}{2}x^2 - x + B $.

The method of undetermined coefficients (I)

Problem:

$ y'' - y = \cos x $

Solution:

Let's first solve $ z'' - z  = 0 $, using the characteristic polynomial method, we know the answer is $ z = Ae^x + Be^{-x} $.

It remains to find just one function that fits $ w'' - w = \cos x $, then we know $ y = w + z $ and have the degree of freedom we need.

At this point, the so called "method of undetermined coefficients" is really just about guessing the solution, let's guess our solution is $ w = \cos x $ and see

$ (\cos x)'' - \cos x = -2 \cos x $.

So now the solution is obviously $ w = \frac{-1}{2}\cos x $

The full solution is then $ Ae^x + Be^{-x} - \frac{1}{2}\cos x $.

Second order constant coefficients (III)

Problem:

2y'' - 3y' + 6y = 0

Solution:

The solutions to the characteristic polynomial are $ \frac{3 \pm \sqrt{39}i}{4} $

Skipping all the derivation we have already done in the last post, the answer is

$ y = Ae^{\frac{3}{4}}\cos\frac{\sqrt{39}}{4}x + Be^{\frac{3}{4}}\sin\frac{\sqrt{39}}{4}x $

Second order constant coefficients (II)

Problem:

$ y'' + y' + y = 0 $

Solution:

The characteristic polynomial is $ x^2 + x + 1 = 0 $, the solution is $ \frac{-1 \pm \sqrt{1^2 - 4(1)(1)}}{2(1)} = \frac{-1 \pm \sqrt{3} i}{2} $

Therefore the answer is

$ y = Ae^{\frac{-1 + \sqrt{3} i}{2}x} + Be^{\frac{-1 - \sqrt{3} i}{2}x} $

Separating out the real and imaginary parts, we get

$ y = Ae^{\frac{-1}{2}x} e^{\frac{\sqrt{3} i}{2}x} + Be^{\frac{-1}{2}x}e^{\frac{- \sqrt{3} i}{2}x} $

Next we can apply the Euler's formula

$ y = Ae^{\frac{-1}{2}x} (\cos\frac{\sqrt{3}}{2}x + i\sin\frac{\sqrt{3}}{2}x) + Be^{\frac{-1}{2}x}(\cos\frac{\sqrt{3}}{2}x - i\sin\frac{\sqrt{3}}{2}x) $

Grouping the real and imaginary gives:

$ y = Ae^{\frac{-1}{2}x} \cos\frac{\sqrt{3}}{2}x  + Be^{\frac{-1}{2}}\cos\frac{\sqrt{3}}{2}x + Ae^{\frac{-1}{2}}\sin(\frac{\sqrt{3}}{2}x)i - Be^{\frac{-1}{2}x}\sin(\frac{\sqrt{3}}{2}x)i $

Factorize, we get:

$ y = (A + B) e^{\frac{-1}{2}x} \cos\frac{\sqrt{3}}{2}x + (A-B)e^{\frac{-1}{2}}\sin(\frac{\sqrt{3}}{2}x)i $

There are two ways we can make the expression real, either we set $ A = B $, that get rid of the imaginary part, or we can set $ A = -B $ be a pure imaginary number, that essentially get rid of the real part and make the imaginary part real. That gives the answer as:

$ y = C e^{\frac{-1}{2}x} \cos\frac{\sqrt{3}}{2}x + D e^{\frac{-1}{2}}\sin\frac{\sqrt{3}}{2}x $

On repeated roots

Following up with the last post, we wanted to understand why $ xe^{-2x} $ is the solution when $ -2 $ is a repeated root.

Let's start with this, $ (p(x) e^{ax})^{(k)} $

We notice this pattern

$ (p(x)e^{ax})' = e^{ax}(ap(x) + p'(x)) $
$ (p(x)e^{ax})'' = e^{ax}(a^2p(x) + 2ap'(x) + p''(x)) $
$ (p(x)e^{ax})''' = e^{ax}(a^3p(x) + 3a^2p'(x) + 3ap''(x) + p'''(x)) $

Feel like binomial expansion? Let's prove this.

Let $ S(n) $ be the statement $ (p(x)e^{ax})^{(n)} = e^{ax}\sum\limits_{i=0}^{n}(\left(\begin{array}{c}n\\i\end{array}\right)a^i p^{(n-i)}(x)) $

$ S(0) $ is obviously true, now let's assume $ S(k) $ is true and

$ (p(x)e^{ax})^{(k)} = e^{ax}\sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k-i)}(x)) $

Now we differentiate both sides once more, we get

$ (p(x)e^{ax})^{(k+1)} = e^{ax} (\sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k-i)}(x)))' + ae^{ax}\sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k-i)}(x)) $

Don't be shied away from the first term, it is just differentiating polynomials, so adding 1 to the differentiation times and that's all. grouping terms, we will get

$ \begin{eqnarray*} (p(x)e^{ax})^{(k+1)} &=& e^{ax} \sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k+1-i)}(x)) + ae^{ax}\sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k-i)}(x)) \\ &=& e^{ax} \sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k+1-i)}(x)) + ae^{ax}\sum\limits_{i=1}^{k+1}(\left(\begin{array}{c}k\\i - 1\end{array}\right)a^{i-1} p^{(k+1-i)}(x)) \\ &=& e^{ax} \sum\limits_{i=0}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k+1-i)}(x)) + e^{ax}\sum\limits_{i=1}^{k+1}(\left(\begin{array}{c}k\\i - 1\end{array}\right)a^i p^{(k+1-i)}(x)) \\ &=& e^{ax} \left(\begin{array}{c}k\\0\end{array}\right)a^0 p^{(k+1-0)}(x) + e^{ax} \sum\limits_{i=1}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k+1-i)}(x)) + e^{ax}\sum\limits_{i=1}^{k}(\left(\begin{array}{c}k\\i - 1\end{array}\right)a^i p^{(k+1-i)}(x)) + e^{ax}\left(\begin{array}{c}k\\k+1 - 1\end{array}\right)a^{k+1} p^{(k+1-(k+1))}(x) \\ &=& e^{ax} a^0 p^{(k+1)}(x) + e^{ax} \sum\limits_{i=1}^{k}(\left(\begin{array}{c}k\\i\end{array}\right)a^i p^{(k+1-i)}(x)) + e^{ax}\sum\limits_{i=1}^{k}(\left(\begin{array}{c}k\\i - 1\end{array}\right)a^i p^{(k+1-i)}(x)) + e^{ax}a^{k+1} p^{(0)}(x) \\ &=& e^{ax} a^0 p^{(k+1)}(x) + e^{ax} \sum\limits_{i=1}^{k}(\left(\left(\begin{array}{c}k\\i\end{array}\right)+\left(\begin{array}{c}k\\i-1\end{array}\right)\right)a^i p^{(k+1-i)}(x)) + e^{ax}a^{k+1} p^{(0)}(x) \\ &=& e^{ax} a^0 p^{(k+1)}(x) + e^{ax} \sum\limits_{i=1}^{k}(\left(\begin{array}{c}k + 1\\i\end{array}\right)a^i p^{(k+1-i)}(x)) + e^{ax}a^{k+1} p^{(0)}(x) \\ &=& e^{ax} \sum\limits_{i=0}^{k+1}(\left(\begin{array}{c}k + 1\\i\end{array}\right)a^i p^{(k+1-i)}(x)) \\ \end{eqnarray*} $

Well, a little more complicated than I wanted it to be, but I proved the result.

Next, we assume the differential equation has constant coefficient like this $ \sum\limits_{j=0}^{n} c_j y^{(j)} = 0 $, and that $ r $ is a root of multiplicity $ k $. In this case, we claim $ x^m e^{rx} $ is a solution for the differential equation for $ 0 \le m < k $. All we need to do is to put it back to the equation and use our result above:

$ \begin{eqnarray*} & & \sum\limits_{j=0}^{n} c_j y^{(j)} \\ &=& \sum\limits_{j=0}^{n} c_j e^{rx}\sum\limits_{i=0}^{j}(\left(\begin{array}{c}j\\i\end{array}\right)r^i p^{(j-i)}(x)) \\ &=& e^{rx} \sum\limits_{j=0}^{n} c_j \sum\limits_{i=0}^{j}(\left(\begin{array}{c}j\\i\end{array}\right)r^i p^{(j-i)}(x)) \\ &=& e^{rx} \sum\limits_{j=0}^{n} c_j \sum\limits_{i=0}^{j}(\left(\begin{array}{c}j\\i\end{array}\right)r^{j-i} p^{(i)}(x)) \\ &=& e^{rx} \sum\limits_{i=0}^{n} \sum\limits_{j=i}^{n} c_j (\left(\begin{array}{c}j\\i\end{array}\right)r^{j-i} p^{(i)}(x)) \\ &=& e^{rx} \sum\limits_{i=0}^{n} \sum\limits_{j=i}^{n} c_j (\frac{j!}{i! (j-i)!}r^{j-i} p^{(i)}(x)) \\ &=& e^{rx} \sum\limits_{i=0}^{n} \frac{p^{(i)}(x)}{i!} \sum\limits_{j=i}^{n} c_j (\frac{j!}{(j-i)!}r^{j-i} ) \\ \end{eqnarray*} $

Note that the inner summation is really just the characteristics polynomial differentiated $ i $ times evaluated at $ r $, and we know $ r $ is a root of multiplicity $ k $, so for all $ i < k $, the term is zero.

For the other terms, we differentiated $ p(x) $ at least $ k $ times, but $ p(x) = x^m $, so the derivative vanishes, and those terms go to 0 as well. Combining the two facts, the differential equation is satisfied, and therefore it is a solution!

Q.E.D.

The magic step of the fourth equal sign probably need some explanation. One can better visualize that if one programs, it is basically a loop transformation

for (j = 0 to n) { for (i = 0 to j } { } }

is transformed to

for (i = 0 to n) { for (j = i to n } { } }

The easiest way to see why this transformation is valid is by generating the (i, j) tuples created by these loops.

Second order constant coefficients (I)

Problem:

$ y'' + 4y' + 4y = 0 $

Solution:

Note that we do not have the independent variable $ x $, in principle we could have reduce to order to 1, but for this one let's try the solution $ y = e^{rx} $. Putting the solution in the problem, we have

$ r^2e^{rx} + 4re^{rx} + 4e^{rx} = 0 $

Solving, get $ r = -2 $ (repeated).

The repeated solution worried me because then we have only one degree of freedom.

Let's try $ y = xe^{rx} $, $ y' = e^{rx} + rxe^{rx} $, $ y'' = re^{rx} + re^{rx} + r^2xe^{rx} $

$ y'' + 4y' + 4y = (re^{rx} + re^{rx} + r^2xe^{rx}) + 4(e^{rx} + rxe^{rx}) + 4(xe^{rx}) $
$  = e^{rx}((r + r + r^2x) + 4(1 + rx) + 4(x)) $
$  = e^{rx}(2r + 4+ (r^2 + 4r + 4)) $

So if we put $ r = -2 $, the equation still work out, so the solution is

$ Ae^{-2x} + Bxe^{-2x} $.

To be honest, I kind of know $ xe^{-2x} $ is an answer when I see the double root. I know it because of my experience (sort of my teacher told me so thing). But wouldn't it be more satisfying if we figured out the general rule and why the general rule work? Let's investigate in the next post.

Reduction of order (II)

Problem:

$ y'' - y'y = 0 $

Solution:

This time we see $ x $ does not appear in the equation, we can let $ z = y' $ and $ y'' = \frac{dz}{dx} = \frac{dz}{dy}\frac{dy}{dx} = z'z $, note that $ z' $ denote $ \frac{dz}{dy} $.

The equation becomes

$ \begin{eqnarray*} y'' - y'y &=& 0 \\ z'z - z y &=& 0 \\ z'z &=& z y \\ z' &=& y \\ z &=& \frac{y^2}{2} + c_1 \\ y' &=& \frac{y^2}{2} + c_1 \\ \frac{2dy}{y^2 + c_2^2} &=& dx \\ 2 \arctan(y/c_2)/c_2 &=& x + c_3 \\ \arctan(y/c_2) &=& \frac{c_2x + c_4}{2} \\ y/c_2 &=& \tan\frac{c_2x + c_4}{2} \\ y &=& c_2 \tan\frac{c_2x + c_4}{2} \\ y &=& 2 c_5 \tan(c_5x + c_6) \\ \end{eqnarray*} $

Check:

$ \begin{eqnarray*} y &=& 2 c_5 \tan(c_5x + c_6) \\ y' &=& 2 c_5 \sec^2(c_5x + c_6) (c_5) \\ &=& 2 c_5^2 \sec^2(c_5x + c_6) \\ y'' &=& 2 c_5^2 (2 \sec(c_5x + c_6))(\sec(c_5x + c_6) \tan(c_5x + c_6)) (c_5) \\ &=& 4 c_5^3 \sec^2(c_5x + c_6) \tan(c_5x + c_6) \end{eqnarray*} $

Reduction of order (I)

Problem:

$ y'' - y' = x $

Solution:

y is missing from the equation, so we can simply let $ z = y' $ and so $ y'' = z' $. The equation is reduced to

$ \frac{dz}{dx} - z = x $

This is a first order linear equation, by inspection, the integrating factor is $ e^{-x} $, we have

$ \frac{d}{dx}(e^{-x} z) = e^{-x} \frac{dz}{dx} - e^{-x} z = e^{-x} (\frac{dz}{dx} - z) = x e^{-x} $.

Therefore $ e^{-x} z = \int {x e^{-x} dx } = - xe^{-x} - e^{-x} + C_1 $.

Then we have $ z = -x - 1 + C_1e^x $

Last but not least, $ y = \int z dx = -\frac{x^2}{2} - x + C_1e^x + C_2 $.

Integrating Factor

Problem:

$ \frac{2y}{x^2}dx + \frac{1}{x}dy = 0 $.

Solution:

Let $ M(x, y) = \frac{2y}{x^2} $, $ N(x, y) = \frac{1}{x} $, now we have

$ \frac{\partial M}{\partial y} = \frac{2}{x^2} \ne \frac{-1}{x^2} = \frac{\partial N}{\partial x} $

But they look similar, so let's assume there exists $ \mu(x) $ such that

$ \frac{\partial \mu M}{\partial y} = \frac{\mu \partial N}{\partial x} $

Expanding, we get

$ \mu \frac{\partial M}{\partial y} = \frac{\partial \mu}{\partial x} N + \mu \frac{\partial N}{\partial x} $

or

$ \frac{1}{N}(\frac{\partial M}{\partial y} - \frac{\partial N}{\partial x})= \frac{1}{\mu}\frac{\partial \mu}{\partial x}$

so

$ \mu = e^{\int \frac{1}{N}(\frac{\partial M}{\partial y} - \frac{\partial N}{\partial x}) dx} $

Putting in the values, we have

$ \begin{eqnarray*} \mu &=& e^{\int \frac{1}{\frac{1}{x}}(\frac{2}{x^2} - \frac{-1}{x^2}) dx} \\ &=& e^{\int \frac{1}{\frac{3}{x}} dx} \\ &=& e^{3 \log x} \\ &=& x^3 \end{eqnarray*} $

The rest is rather easy now, we need to find $ f(x, y) $ such that $ \frac{\partial f}{\partial x} = x^3 \frac{2y}{x^2} = 2xy $ and $ \frac{\partial f}{\partial y} = x^3 \frac{1}{x} = x^2 $. We found $ f(x, y) = x^2y = c $, therefore the solution to the differential equation is $ y = cx^{-2} $.

Homogenous Equation

Problem:

$ (y^2 - x^2)dx + xydy = 0 $

Solution:

$ \begin{eqnarray*} (y^2 - x^2)dx + xydy &=& 0 \\ (y^2 - x^2)dx &=& -xydy \\ \frac{dy}{dx} &=& \frac{x^2 - y^2}{xy} \\ &=& \frac{x}{y} - \frac{y}{x} \\ \end{eqnarray*} $

Now let $ z = \frac{y}{x}, \frac{dz}{dx} = \frac{d}{dx}(y)\frac{1}{x} + y\frac{d}{dx}\frac{1}{x} = \frac{dy}{dx}\frac{1}{x} - \frac{y}{x^2} $, which implies $ \frac{dy}{dx} = x\frac{dz}{dx} + z $. Putting these back, we have

$ \begin{eqnarray*} \frac{dy}{dx} &=& \frac{x}{y} - \frac{y}{x} \\ x\frac{dz}{dx} + z &=& \frac{1}{z} - z \\ x\frac{dz}{dx} &=& \frac{1}{z} - 2z \\ \frac{zdz}{1-2z^2} &=& \frac{dx}{x} \\ \frac{-1}{4}\log(1-2z^2) &=& \log(x) + C_1 \\ \log(1-2z^2) &=& -4\log(x) + C_2 \\ 1-2z^2 &=& C_3x^{-4} \\ 1-2(\frac{y}{x})^2 &=& C_3x^{-4} \\ x^2-2y^2 &=& C_3x^{-2} \\ y^2 &=& \frac{1}{2}(x^2 - Cx^{-2}) \\ \end{eqnarray*} $

For a quick check, we have:

$ \begin{eqnarray*} y^2 &=& \frac{1}{2}(x^2 - Cx^{-2}) \\ 2y \frac{dy}{dx} &=& x + Cx^{-3} \\ &=& x + \frac{Cx^{-2}}{x} \\ &=& x + \frac{x^2 - 2y^2}{x} \\ &=& 2x - \frac{2y^2}{x} \\ \frac{dy}{dx} &=& \frac{x}{y} - \frac{y}{x} \\ \end{eqnarray*} $

Orthogonal Trajectories

Problem:

Find the orthogonal trajectories to the curves $ y = x^2 + c $.

Solution:

First, we express the parabolas as differential equations

$ \frac{dy}{dx} = 2x $

Next, the slope of the tangent lines for the orthogonal trajectories is determined, so we solve this

$ \frac{dy}{dx} = \frac{-1}{2x} $

This is separable, and therefore the solution is simply

$ y = \frac{-1}{2}sgn(x)\ln|x| + c $

Exact differential equation

Problem:

$ 3x^2y dx + x^3dy = 0 $

Solution:

We have $ \frac{\partial 3x^2y}{\partial y} = 3x^2 = \frac{\partial x^3}{\partial x} $, therefore, we can find $ f(x, y) $ such that $ \frac{\partial f}{\partial x} = 3x^2y $ and $ \frac{\partial f}{\partial y} = x^3 $

It is quite obvious that $ f(x, y) = x^3y + c $.

The answer to the differential equation is therefore $ y = cx^{-3} $.

As a quick check, we have $ \frac{dy}{dx} = -3cx^{-4} $, therefore

$ 3x^2y dx + x^3dy = 3x^2(cx^{-3})dx + x^3(-3cx^{-4})dx = 0 $.

First order linear differential equation

Problem:

$ y' + 2xy = x $

Solution:

This is a first order linear differential equation, let's consider the derivative of the following form:

$ (p(x)e^{q(x)})' = p'(x)e^{q(x)} + p(x)q'(x)e^{q(x)} = e^{q(x)}(p'(x) + p(x)q'(x)) $

Now we set $ p(x) = y $, $ q'(x) = 2x $, we can multiply both side by $ e^{x^2} $ to get

$ \begin{eqnarray*} e^{x^2}(y' + 2xy) &=& xe^{x^2} \\ (e^{x^2}y)' &=& xe^{x^2} \\ (e^{x^2}y) &=& \int{xe^{x^2}dx} \\ &=& \frac{1}{2}\int{e^{x^2}d(x^2)} \\ &=& \frac{1}{2}e^{x^2} + C \\ y &=& \frac{1}{2} + Ce^{x^{-2}} \\ \end{eqnarray*} $

Coupled differential equation

Today morning, I am having fun solving this coupled differential equation pair.

$ \dot{X} = \left[\begin{array}{cc} 1 & -4 \\ 4 & 1\end{array}\right] X + \left[\begin{array}{c} 4t + 9e^{6t}\\-t+e^{6t} \end{array}\right] $

With initial conditions $ X(0) = \frac{1}{17}\left[\begin{array}{c}21\\35\end{array}\right] $

The key to solving this puzzle is to decouple the equations, to do that, we diagonalize the matrix. I get this handy equation:

$ \left[\begin{array}{cc} 1 & -4 \\ 4 & 1\end{array}\right] = \frac{1}{2}\left[\begin{array}{cc} 1 & 1 \\ -i & i\end{array}\right]\left[\begin{array}{cc} 1 +4i & 0 \\ 0 & 1 - 4i\end{array}\right]\left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right] $

We substitute this back to the equation, we see something interesting.


$ \begin{eqnarray} \dot{X} &=& \frac{1}{2}\left[\begin{array}{cc} 1 & 1 \\ -i & i\end{array}\right]\left[\begin{array}{cc} 1 +4i & 0 \\ 0 & 1 - 4i\end{array}\right]\left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right] X + \left[\begin{array}{c} 4t + 9e^{6t}\\-t+e^{6t} \end{array}\right] \\ \left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]\dot{X} &=& \frac{1}{2}\left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]\left[\begin{array}{cc} 1 & 1 \\ -i & i\end{array}\right]\left[\begin{array}{cc} 1 +4i & 0 \\ 0 & 1 - 4i\end{array}\right]\left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right] X + \left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]\left[\begin{array}{c} 4t + 9e^{6t}\\-t+e^{6t} \end{array}\right] \\ \left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]\dot{X} &=& \left[\begin{array}{cc} 1 +4i & 0 \\ 0 & 1 - 4i\end{array}\right]\left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right] X + \left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]\left[\begin{array}{c} 4t + 9e^{6t}\\-t+e^{6t} \end{array}\right] \\ \end{eqnarray} $

It is now obvious that we should let $ Y = \left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]X $. That would then give

$ \begin{eqnarray} \dot{Y} &=& \left[\begin{array}{cc} 1 +4i & 0 \\ 0 & 1 - 4i\end{array}\right]Y + \left[\begin{array}{cc} 1 & i \\ 1 & -i\end{array}\right]\left[\begin{array}{c} 4t + 9e^{6t}\\-t+e^{6t} \end{array}\right] \\ \end{eqnarray} $

Now we have decoupled the equations. Let $ Y = \left[\begin{array}{c}y_1\\y_2\end{array}\right] $, we have two separate equations instead of a coupled pair, which can be solved separately.

$ \begin{eqnarray} y_1'(t) &=& (1 + 4i)y_1(t) + (4t + 9e^{6t}) + i(-t + e^{6t}) \\ y_2'(t) &=& (1 - 4i)y_2(t) + (4t + 9e^{6t}) - i(-t + e^{6t}) \end{eqnarray} $

Just so we are not lost, we also let $ X = \left[\begin{array}{c}x_1\\x_2\end{array}\right] $, and therefore we have these relations:

$ \begin{eqnarray} y_1 &=& x_1 + ix_2 \\ y_2 &=& x_1 - ix_2 \\ \end{eqnarray} $

Let's try to solve it with Laplace's transform. For the first equation, taking Laplace's transform on both sides give this:

$ \begin{eqnarray} y_1'(t) &=& (1 + 4i)y_1(t) + (4t + 9e^{6t}) + i(-t + e^{6t}) \\ &=& (1 + 4i)y_1(t) + (4-i)t + (9+i)e^{6t} \\ \mathcal{L}(y_1'(t)) &=& \mathcal{L}((1 + 4i)y_1(t) + (4-i)t + (9+i)e^{6t}) \\ sY_1(s) - y_1(0) &=& (1 + 4i)Y_1(s) + \frac{4-i}{s^2} + \frac{9+i}{s - 6} \\ sY_1(s) - (1 + 4i)Y_1(s) &=& y_1(0) + \frac{4-i}{s^2} + \frac{9+i}{s - 6} \\ (s - 1 - 4i)Y_1(s) &=& \frac{21 + 35i}{17} + \frac{4-i}{s^2} + \frac{9+i}{s - 6} \\ Y_1(s) &=& \frac{21 + 35i}{17(s - 1 - 4i)} + \frac{4-i}{s^2(s - 1 - 4i)} + \frac{9+i}{(s - 6)(s - 1 - 4i)} \\ \end{eqnarray} $

Without repeating, we get, similarly, that

$ \begin{eqnarray} Y_2(s) &=& \frac{21 - 35i}{17(s - 1 + 4i)} + \frac{4+i}{s^2(s - 1 + 4i)} + \frac{9-i}{(s - 6)(s - 1 + 4i)} \\ \end{eqnarray} $

The rest is just ugly partial fractions and inverse Laplace transform.