Day 14

Tuesday, March 28, 2017

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

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*)*u*_{t}(*x*, *t*) = (*A*(*x*)*K*_{0}(*x*)*u*_{x}(*x*, *t*))_{x}.

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

*c*(*x*)*ρ*(*x*)*A*(*x*)*φ*(*x*)*h*^{′}(*t*) = (*A*(*x*)*K*_{0}(*x*)*φ*^{′}(*x*))^{′}*h*(*t*)

which separates as

(*h*^{′}(*t*))/(*h*(*t*)) = ((*A*(*x*)*K*_{0}(*x*)*φ*^{′}(*x*))^{′})/(*c*(*x*)*ρ*(*x*)*A*(*x*)*φ*(*x*)) = − *λ*

giving

*h*^{′}(*t*) = − *λ**h*(*t*),

(*A*(*x*)*K*_{0}(*x*)*φ*^{′}(*x*))^{′} + *λ**c*(*x*)*ρ*(*x*)*A*(*x*)*φ*(*x*) = 0.

Let us assume the bar has uniform material composition, so that *K*_{0}
*ρ* and *c* do not depend on location - say *K*_{0}(*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.

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

(Python code to draw bar here)

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!

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?

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.

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:

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".

5 very definitive things

in a slightly generalized sense

up to multiplication by a constant

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:

∫^{L}_{0}*φ*_{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*))/(√(∫^{L}_{0}*φ*(*x*)^{2}*σ*(*x*)*dx*)).

Then if the "scalar product"

∫^{L}_{0}*φ*_{1}(*x*)*φ*_{2}(*x*)*σ*(*x*)*dx*

is small *compared to 1*, *φ*_{1} and *φ*_{2} can be
considered approximately orthogonal.