Day 14

Tuesday, March 28, 2017

# Ch 5: Sturm-Liouville problems, cont'd

Let's generalize the problems we've been looking at so far.

## Heat flow in a bar of non-uniform cross-section

Going way back to Ch. 1, we have the following PDE for the temperature profile u(x, t) in a slender bar, allowing all the material properties all to be non-constant:

c(x)ρ(x)A(x)ut(x, t) = (A(x)K0(x)ux(x, t))x.

Setting u(x, t) = φ(x)h(t) this says

c(x)ρ(x)A(x)φ(x)h(t) = (A(x)K0(x)φ(x))h(t)

which separates as

(h(t))/(h(t)) = ((A(x)K0(x)φ(x)))/(c(x)ρ(x)A(x)φ(x)) =  − λ

giving

h(t) =  − λh(t),

(A(x)K0(x)φ(x)) + λc(x)ρ(x)A(x)φ(x) = 0.

Let us assume the bar has uniform material composition, so that K0 ρ and c do not depend on location - say K0(x) = 1, ρ(x) = 1, c(x) = 1, but we allow the cross-sectional area, A(x), to vary.

Also, let's assume we will hold the left-hand end (x=0) at 0 degrees, and insulate the opposite end (x=L) so the flux is zero at x=L.

Then the spatial eigenproblem becomes:

(A(x)φ(x)) + λA(x)φ(x) = 0,

φ(0) = 0,  φ(L) = 0.

If we can find some solutions of this ODE BVP, then they are spatial eigenfunctions of the PDE. A priori they might not be orthogonal, though, and it might be possible that we'll have only a limited number of them. (The Sturm-Liouville theory will answer these worries.)

Let us pick some specific varying cross-section function, A(x) and then use an ODE solver in Python to try and find the eigenvalues and eigenfunctions for the corresponding ODE for φ(x) using the "shooting" method I will describe.

### Exercise 1

Suppose we choose A(x) = 1.1 + cos(2πx)/(L)

What do you think the slowest-decaying (i.e. smallest λ) eigenfunction will look like in this problem? Think about it, then draw your guess on the board, and sign it!

### Exercise 2

Will the λ1 for this non-uniform bar be larger (i.e. mode decays faster) or smaller (i.e. mode decays slower) than the λ1 for the uniform bar?

### Numerically compute the first eigenfunction and the corresponding eigenvalue

Numerical ODE solvers are all set up to solve systems of 1st order ODES. So we'll need to convert our 2nd order ODE to that form.

First expand the first term of the ODE using the product rule as:

(Aφ) = Aφ′′ + Aφ.

Then continue by defining φ′ = ψ, say.

GOT TO HERE Day 13.

Specifically, let us find the first (smallest) several eigenvalues and make a plot the corresponding eigenfunctions.

Here's my code to do that

## Shooting to find eigenvalues and eigenfunctions

### Smallest eigenvalue

Our first shots:

ringland@blue:~/public_html/418$py non-uniform_bar.py lambda = 0.0 : phi'(L) = 1.000000444587765 lambda = 0.01 : phi'(L) = -2.5218656367781027 We seem to have captured the smallest eigenvalue somewhere between 0 and 0.01. Let's look closer ... ringland@blue:~/public_html/418$ py non-uniform_bar.py
lambda = 0.0                 : phi'(L) = 1.000000444587765
lambda = 0.00025             : phi'(L) = 0.8808779967983755
lambda = 0.0005              : phi'(L) = 0.7635397347055736
lambda = 0.00075             : phi'(L) = 0.6479696072306252
lambda = 0.001               : phi'(L) = 0.5341519836659785
lambda = 0.00125             : phi'(L) = 0.42207088354630423
lambda = 0.0015              : phi'(L) = 0.31171066959066557
lambda = 0.00175             : phi'(L) = 0.20305583065389565
lambda = 0.002               : phi'(L) = 0.09609071981800561
lambda = 0.00225             : phi'(L) = -0.0091995598926174
lambda = 0.0025              : phi'(L) = -0.1128306966423111
lambda = 0.00275             : phi'(L) = -0.21481770101420766
[remainder omitted]

So our smallest eigenvalue looks to be λ1 ≈ 0.00225. Therefore this is the picture to compare your prediction sketch of Day 13:

### How about a second eigenvalue?

We'll try going up to lambda = 0.08:

ringland@blue:~/public_html/418$py non-uniform_bar.py lambda = 0.0 : phi'(L) = 1.000000444587765 lambda = 0.01 : phi'(L) = -2.5218656367781027 lambda = 0.02 : phi'(L) = -4.051481202064394 lambda = 0.03 : phi'(L) = -4.279579716989934 lambda = 0.04 : phi'(L) = -3.7201642805170243 lambda = 0.05 : phi'(L) = -2.746567965301901 lambda = 0.06 : phi'(L) = -1.6213808523489763 lambda = 0.07 : phi'(L) = -0.521142817358817 lambda = 0.08 : phi'(L) = 0.44341238297076263 Good! We have passed through zero-slope at x=L again. Zooming in a bit ... ringland@blue:~/public_html/418$ py non-uniform_bar.py
lambda = 0.07                : phi'(L) = -0.521142817358817
lambda = 0.07200000000000001 : phi'(L) = -0.31479480787128866
lambda = 0.07400000000000001 : phi'(L) = -0.11465232772346716
lambda = 0.076               : phi'(L) = 0.07873781172935942
lambda = 0.078               : phi'(L) = 0.2648969971065854
lambda = 0.08                : phi'(L) = 0.44341238297076263

Looks like λ2 ≈ 0.0752:

We can continue to increase lambda to find additional eigenvalues. Here 50 values between 0.0753 and 1:

Looks like we passed over several eigenvalues in that range. And we can go higher still. The next set of shots shows there is an eigenvalue λ ≈ 51.12:

ringland@blue:~/public_html/418\$ py non-uniform_bar.py
lambda = 51.0                : phi'(L) = -0.16168482969249998
lambda = 51.02               : phi'(L) = -0.1339811444223814
lambda = 51.04               : phi'(L) = -0.10617754016583797
lambda = 51.06               : phi'(L) = -0.078296680951926
lambda = 51.08               : phi'(L) = -0.05035974923012558
lambda = 51.1                : phi'(L) = -0.022388623013249467
lambda = 51.12               : phi'(L) = 0.00559413827806367
lambda = 51.14               : phi'(L) = 0.033568092967984964
lambda = 51.16               : phi'(L) = 0.06150911907898296
lambda = 51.18               : phi'(L) = 0.08939673734343624
lambda = 51.2                : phi'(L) = 0.11720962485962341

It would be nice to compare and contrast with the eigenthings for a uniform bar, and see if our eigenvalues and eigenfunctions (especially λ1 and φ1(x)) for this particular non-uniform bar "make sense".

# Things we can say about all Sturm-Liouville problems

5 very definitive things

# Proofs of three of the things

## Eigenfunctions are orthogonal

in a slightly generalized sense

## There is only one eigenfunction for each eigenvalue

up to multiplication by a constant

# Back to our non-uniform bar

## Can we verify that the eigenfunctions mutually orthogonal?

Can we also check to see if your eigenfunctions are orthogonal? (Just one pair, let's say.)

We will see that we actually have to modify the definition of orthogonality slightly to make this work:

L0φ1(x)φ2(x)σ(x)dx = 0

where σ(x) is a weighting function specific to the DE we're solving.

If we have two eigenfunctions and the weighting function all evaluated at a grid of x-values (the same grid for all), then the integral of their product can be approximated by forming the elementwise product of the three vectors and summing.

Note on checking orthogonality: We are computing approximations to the eigenfunctions, so they cannot be expected to be exactly orthogonal. So when we compute the (weighted) integral of their product, which you expect to be zero, you have to have some sense of how small is close enough to zero. Whatever this measure is, it should account for the fact that eigenfunctions are determined only up to scale. Thus your measure of their (lack of) orthogonality should be invariant (unchanged) if either of the eigenfunctions is rescaled. One way of doing this is to normalize all your eigenfunctions (analogous to making them "unit vectors") as follows:

φ(x) → (φ(x))/((L0φ(x)2σ(x)dx)).

Then if the "scalar product"

L0φ1(x)φ2(x)σ(x)dx

is small compared to 1, φ1 and φ2 can be considered approximately orthogonal.