Methods in Applied Mathematics

(1)

Methods in Applied Mathematics

I-Liang Chern

Department of Mathematics National Taiwan University

(2)

(3)

Dimensional Analysis

1.1 The program of applied mathematics

The program of applied mathematics consists of the following procedures.

• Observe a phenomenon, determine variables (quantities) that we are concerned with. Phys-ical quantities have dimension such as “length”L, “time”T, “mass”M, “charge”e. Quantify the variables by choosing “UNIT”. Many physical quantities are derivatives of certain funda-mental dimensions. For instance, the dimension of energy is M L2T−2. The dimension of a variable is denoted by [ · ].

• Modeling: Find relations among quantities 1. determine dependent / independent variables 2. determine parameters

3. The relation can be found based on basic physical laws and/or dimensional analysis 4. Relations are in the forms of algebraic relations, ordinary/partial differential equations,

integral equations, ...

5. The key is “In a relation, each term has to be of the same dimension.” • Simplification and Reduction

1. Reduce numbers of parameters by dimensional analysis

2. Reduce numbers of unknowns (dependent variables) by perturbation methods 3. Sensitivity analysis to find important quantities

• Find solutions by analytic, computational methods

• Perform experiments to check solutions, modify models if needed. 3

(7)

4 CHAPTER 1. DIMENSIONAL ANALYSIS

1.2 Dimensional Analysis

Population model dP dt = rP − bP 2_.

Here, r is the growth rate, b is the competition rate per each individual. It means that the species compete for the same resource and causes the reduction of population. Thus, this is an environmen-tal constraint. The initial population is P0. In this model, the quanities are P , t, P0, r and b. t is the

independent variable, P the depedent variable, P0, r and b are parameters. The dimensions of them

are [t] = T, [P ] = P, [P0] = P, [r] = 1 T, [b] = 1 T P.

We see that there are only two fundamental dimensions T and P . We can rescale this equation by a proper scaling as the follows. We write the equation by

dP d(rt) = P − b rP 2 _{= P} 1 − P K , K = r b. Thus, we can introduce the following dimensionless variables

t0 = rt, P0 = P K, P 0 0 = P0 K, then the equation is equivalent

dP0 dt0 = P

0

(1 − P0), P0(0) = P₀0.

This equivalent equation depends only on 3 dimensionless quantities P0, t0 and P₀0, instead of 5 quantities.

A falling object 1. Quantities

(a) independent variables: t. (b) dependent variables: x (c) parameters: g, V , m, h 2. Relations Newton’s law:

m¨x = −mg, x(0) = h, ˙x(0) = V, 3. Dimensional Analysis

(8)

1.2. DIMENSIONAL ANALYSIS 5 (b) [t] = T ,

(c) [g] = L/T2. (d) [V ] = L/T

We will find a relation between x, t, g, h, V . Notice that there are only two fundamental dimensions, namely, L and T .

4. Scaling: A natural way is to introduce h and V as our fundamental dimensions, which are the observed (or controlled) parameters. Then we rescale x0 = x/h, v0 = v/V , t0 = t/T = t/(h/V ). With these, the equation becomes

d2x0

dt02 = −gh −1

T2 = −ghV−2= − 1 Fr2. The initial conditions are

x0(0) = 1,dx

0

dt0(0) = 1.

Here, we introduce a dimensionless parameter Fr = √V

gh =

inertial force gravitational force,

called Froude number, which also appears in water wave theory. Then the dimensionless variables are (x0, t0) and the Froude number Fr. Thus, the number of variables and parameters are reduced.

The solution of the above dimensionless equation is x0(t0) = 1 + t0−1

2 t02 Fr2. The critical time that the particle starts to fall is

dx0 dt0(t 0 ) = 0. That is, t0 = Fr2.

Thus, the larger the Froude number is (i.e. inertia force >> gravitation force), the longer time the particle starts to fall. Notice that if we drop an object from height h. The speed it touches the ground is√2gh.

Notice that h can be the length we are interested in, needs not be the initial height. Further, the dimensionless parameters are not uniquely chosen. For instance, another way is to choose v0 = V T /h = V /√gh. Then the rescaled system becomes

d2_x0

(9)

6 CHAPTER 1. DIMENSIONAL ANALYSIS with the initial conditions

x0(0) = 1, dx

0

dt0(0) = v 0_,

a relation involve the three dimensionless quantities x0, t0, v0.

People usually use the first scaling because (1) the initial velocity V is something we can con-trol, it is natural to scale velocity in terms of V ; (2)the parameter Fr appears in the equation, we can study solution property (bifurcation, for instance) in terms of the parameter without solving it directly, whereas the solution property in terms of initial condition usually requires solving the equation first.

Body mass index (BMI) The BMI is defined to be BMI = w

h2,

where w is the weight in kilogram and h is the height in meter. What is the dimension of BMI? Why it measures the amount of fat of your body?

Rain droplet Why the viscous force of an object is propotional with the its area? With a fixed density, why the body force of the object is propotional to its volume? Can you answer what size of the raindrop will fall?

Projectile 1. Quantities:

• independent variable t • Dependent variables: z,

• Parameters: m, g, R (earth radius), V initial velocity 2. The relation: md 2_z dt2 = −G M m (z + R)2 This implies d2z dt2 = − R2g (R + z)2 where g = GM_R2 . z(0) = 0,dz dt(0) = V.

(10)

1.2. DIMENSIONAL ANALYSIS 7 3. Dimensional Analysis: the dimensions of the quantities are

[t] = T, [V ] = LT−1, [z] = L, [g] = LT−2, [R] = L.

There are 5 quantities, but only 2 fundamental dimensions involved. Thus, we expect that there are 3 dimensionless quantities. Let us suppose the dimensionless quantity π is a combi-nation of the 5 quantities in the following form:

[π] = [tα1_zα2_Rα3_Vα4_gα5_]

= Tα1−α4−2α5_Lα2+α3+α4+α5

In order to have π to be dimensionless, i.e. [π] = 1, the admissible exponents leading to the following equations:

α1− α4− 2α5 = 0

α2+ α3+ α4+ α5 = 0.

Let us use α3, α4 and α5as free parameters. We then obtain three sets of solutions:

• α₃ = −1, α4= α5 = 0. This gives α2 = 1 and α1 = 0. We then get π1= z/R.

• α₄ = 1, α3 = α5 = 0. This gives α1 = −α2 = α4= 1. We get π2= tz−1V .

• α5 = 1, α3 = α4 = 0. This gives α1 = 2, α2 = −1. We get π3= t2z−1g.

Alternatively, we can use

πx = z/R, πt= tV /R, π =

gR V2.

Then the differential equation can be written as d2πx dπ2 t = R V2 d2z dt2 = − gR V2 1 (1 + πx)2 = − π (1 + πx)2 .

with πx(0) = 0 and dπ_dπx_t(0) = 1. This is the dimensionless relation, which is a relation among

πx, πtand π. So, numbers of variables are reduced from 5 to 3.

The general strategy to select dimensionless variables are • rescale the spatial variable z by z/R

• rescale the temporal variable t by t/(R/V ).

(11)

8 CHAPTER 1. DIMENSIONAL ANALYSIS Diffusion process Consider an explosion process in the space. We are interested in the heat propagation of this exploration. The variables and parameters of this process are:

• independent variables: x and t;

• dependent variables: u (temperature), q heat flux;

• parameters: energy released at the exploration center, k (conductivity),c (heat capacity at constant volume).

Relations: based on conservation of energy, we have d dt Z Ω cu dx = Z ∂Ω q · (−n) dS This yields Z Ω cutdx = − Z Ω ∇ · q dx This is valid for any arbitrary Ω. Thus, we get

cut+ ∇ · q = 0.

The heat flux is related to the gradient of the temperature. This is the Fourier law: q = −k∇u.

Plug this into trhe above equation, The heat equation becomes ut= ∇ ·

k

c∇u = D 4 u.

Here, the notation 4 = ∇2 is the Laplacian. The dimensions of these quantities are [u] = Θ, [c] = EΘ−1L−3, [q] = ET−1L−2, [k] = ET−1L−1Θ−1, [D] = L2T−1.

In the diffusion equation, the only quantities are x, t, u and D. There are three fundamental dimensions: L, T and Θ. Thus, there is only one dimensionless quantity left. Suppose it is π = tα1_xα2_uα3_Dα4_{. In order to be dimensionless, we get}

π1 = x √ Dt Another one is π2 = uc e x 3 ₌ uc e π 3 1(Dt)3/2

Suppose we have a relation: π2 = g(π1). This leads to a relation which can express the unknown u

in terms of the independent variable x, t and the parameter D: u = e

c(Dt)

−3/2

g(√x

(12)

1.2. DIMENSIONAL ANALYSIS 9 Using dimensional analysis, we can find possible relation such as (1.1), which says that (1) the diffusion length scale x is√Dt; (2) the maximum of the heat decays like (Dt)−3/2in three space dimensions. We obtain this without solving the PDE at all, simply a guess from dimensional analy-sis.

There is another technique that used commonly. We can plug this ansatz into equation. Then we will get an ODE for g. Eventually we will get g(π1) = e−π

2 1_.

Fluid Mechanics The incompressible flows are governed by the following Navier-Stokes equa-tions:

∇ · u = 0,

ρ(ut+ u · ∇u) + ∇p = µ 4 u.

Here, u is the velocity, ρ the density, p the pressure, µ the viscosity. We consider the domain has length scale L and time scale T , or equivalent the velocity has scale U = L/T . Notice that what we usually measures are the typical spatial scale L and the velocity U . The time scale T is a derived dimension. The dimension of [ρ] = M/L−3. It is equivalent to use [ρ] or M . It is more natural to use [ρ] for fluids. The dimension of pressure is

[p] = M U T L2 =

ρL3U

(L/U )L2 = ρU 2_.

By matching the dimensions of the Navier-Stokes equation, we see that the dimension of the vis-cosity µ should be

[µ] = ρU2L−1U−1L2 = ρU L.

Now, we rescale x = Lx0, u = U u0. Hence the time scale T = L/U and time is rescaled by t = t0T = t0L/U . The pressure p is rescaled as p = p0ρU2.

With these, the momentum equation becomes ρU2 L u 0 t0 + u0· ∇0u0+ ∇0p0 = µ U L2 4 0_u0_. Or u0_t0 + u0· ∇0u0+ ∇0p0 = 1 Re4 0_u0

Here, the dimensionless quantity

Re = inertia force viscous force = ρU2L−1 µU L−2 = ρU L µ = U L ν ,

called Renold’s number. The independent variables are (x0, t0), the dependent variables are (u0, p0). The parameter is Re. The parameters (ρ, L, U ) are scaled out. The advantages of these scaling is that the typical length and vlocity are scaled out. Thus, we can perform experiments for small size vehicles at low speed, yet we can get the same results for large size vehicles at high speed.

(13)

10 CHAPTER 1. DIMENSIONAL ANALYSIS

1.3 The Buckingham Pi Theorem

Suppose we have a physical relation

f (q1, · · · , qm) = 0.

A physical relation is called unit free if it is independent of the particular set of units we choose for q1, · · · , qm. “Time” is a dimension, while “year”, “month”, “second” are units of the dimension

“time. There are some fundamental dimensions such as length, time, mass, charge. Suppose there are n fundamental dimensions L1, · · · , Ln, n < m involved in the quantities q1, · · · , qm. Each

qnatity qi has dimension

[qi] = La₁1i· · · Lanni

The matrix (aij)n×mis called the dimension matrix. Suppose the matrix (aij) has rank r. Then we

can determine m − r dimensionless quantities πk= Y i qαik i , k = 1, · · · , m − r, and a relateion g(π1, · · · , πm−r) = 0

such that f = 0 if and only if g = 0. This is called the Buckingham’s π theorem. Proof. Let us first look for a dimensionless quantity π:

π =Y i qαi i . The dimension of π is [π] =Y i   Y j Laji j   αi =Y j L P iajiαi j

In order to have π to be dimensionless , we should require

m

X

i=1

ajiαi = 0, for all j = 1, · · · , n

Since the rank of the dimension matrix A is r, the kernel space of A has dimension m − r. Indeed, let A = V ΣUT be the singular value decomposition of A, where U ,V are orthogonal matrices and Σ has the form

Σ =

(Σ1)r×r 0r×(m−r)

0(n−r)×r 0(n−r)×(m−r)

where Σ1 is a diagonal and invertible matrix. Suppose U = (α1, · · · , αm), then the null space of

(14)

1.3. THE BUCKINGHAM PI THEOREM 11 Q

iq αk,i

i , k = 1, · · · , m. Here, αk= (αk,1, · · · , αk,n)T. From Aαk= 0 for k = r + 1, · · · , m, we

get πkwith k = r + 1, · · · , m are dimensionless.

To find a relation in terms of the dimensionless variables πr+1, · · · , πm, which is also equivalent

to f = 0, we first notice that there is an 1-1 corresponding between (q1, · · · , qm) and (π1, · · · , πm)

because U is non-singular. Thus, for any (π1, · · · , πm), we define

g(π1, · · · , πm) = f (q1, · · · , qm).

To show that g depends only on (πr+1, · · · , πm), we will use the unit free property of f = 0. Let

us rescale the fundamental dimensions by ¯ Lj = λjLj, j = 1, · · · , n. Then qiis scaled to ¯ qi = Y j λaji j qi

The property unit free of f = 0 means that

f (¯q1, · · · , ¯qm) = 0 ⇔ f (q1, · · · , qm) = 0 for any λj > 0., j = 1, ..., n

As translated to g, we get

f (¯q1, · · · , ¯qm) = 0 ⇔ f (q1, · · · , qm) = 0

⇔ g(π1, · · · , πm) = 0 ⇔ g(¯π1, · · · , ¯πm) = 0.

Because πr+1, · · · , πmare dimensionless, we have ¯πr+1 = πr+1, · · · , ¯πm= πmfor any choices of

λj > 0, j = 1, ..., n. We will choose λj such that the dimensions

¯

π1 = · · · = ¯πr = 1.

To see these λj can be found, we derive the conditions for λj as the follows.

¯ πk = Y i ¯ qαk,i i = Y i Y j λajiαk,i j Y i qαk,i i ! = Y j λ P iajiαk,i j Y i q_iαk,i ! = exp   X j (X i ajiαk,i)µj  πk

(15)

12 CHAPTER 1. DIMENSIONAL ANALYSIS Recall that AU = V Σ, let us write Aαk = σkβk, k = 1, · · · , r, σk 6= 0. Then taking log of the

above equations, we get

0 =

n

X

j=1

βk,jµj + ln πk, k = 1, · · · , r.

This set of equations is always solvable because (β1, · · · , βr) are orthogonal. Once we determine

µ, hence λ, then

g(¯π1, · · · , ¯πm) = g(1, · · · , 1, ¯πr+1, · · · , ¯πm).

1.4 Applications of dimensional analysis

1.4.1 Use dimensional analysis for finding possible relation

Gauss-Bonnet Theorem In geometry, the principal curvature κi has dimension L−1. For two

dimensional surfaces, the Gaussian curvature K = κ1κ2 has dimension L−2. Thus, consider a

closed surface Σ, the quantity

Z

Σ

K dS

is dimensionless, because the surface element dS has dimension L2. Thus,R

ΣK dS is scale

inde-pendent. A more ambitious guess is that this quantity is independent of a continuous variation of Σ. That is,R_ΣK dS is a topological quantity. We can then compute this quantity for sphere, torus, etc and get the following conjecture:

Z

Σ

K dS = 2πχ(Σ),

where χ(Σ) is the Euler characteristic of Σ. For orientable compact closed surfaces, χ(Σ) = 2−2g, where g is the genus of the surface. Any such a surface is topologically equivalent to a sphere with some handles attached, and g is the number of the handles.

1.4.2 *Other Applications of Dimensional Analysis in the estimates of PDEs

Poincar`e inequality In PDEs, we usually need to estimate the Lp norm of our variables. The estimates are usually in the forms of inequalities, which should be scale invariant. To find such inequalities, an important technique is the dimensional analysis. We define the dimension of kukLpk

in a small region in Rnby1

[kukLp] = (UpLn)1/p= U Ln/p

The dimension of kDukLpis

[k∇ukLp] = U L p Ln 1/p = U Ln/p−1. 1

(16)

1.4. APPLICATIONS OF DIMENSIONAL ANALYSIS 13 An important inequality is the Poincar`e’s inequality. It relates kuk_L2 and k∇uk_L2 in a bounded

domain with zero boundary condition:

kuk_L2 ≤ Ck∇uk_L2

By comparing the dimension, we see that

U Ln/2= [C]U Ln/2−1

This shows that [C] = L, which is indeed the diameter of the domain.

Sobolev’s inequality We want to know the relation between kuk_Lp∗ and kDukLp.

Theorem 1.1. Assume 1 ≤ p < n, there exists a constant C such that kuk_Lp∗_(Ω)≤ CkDuk_Lp_(Ω)

wherep∗ = np/(n − p). Remark.

• From dimension analysis, we have

[kuk_Lp∗] = U Ln/p ∗

[kDukLp] = U L−1+n/p

The inequality can hold only when the dimensions on both sides match. This leads to n

p∗ =

n p − 1. This gives the formula p∗.

• The above formula says that u can gain some integrability (i.e. p∗ > p) if its derivative Du is already p-integrable. From [Du] = [u] − 1, we see that if Du ∈ Lp, which means that [Du]Ln/p = [u]L−1+n/pis finite as L → 0. Thus, there is some more room of integrability of u.

• If p > n, the constraint is even stronger, we should have p = ∞. In fact, we have kuk_C0,1−n/p ≤ CkDukLp

for all u ∈ C_c∞.

(17)

14 CHAPTER 1. DIMENSIONAL ANALYSIS Interpolation Formula In nonlinear PDEs, it is commonly to estimate the norm of a medium derivatives in terms of high derivatives and low derivatives. This is the interpolation formula. Theorem 1.2. Let 0 ≤ j < m and 1 ≤ q, r ≤ ∞. Then for smooth function u in a bounded domain in Rn, we have kDjukLp ≤ CkDmukα_Lrkuk1−α_Lq , where 1 p = j n+ α 1 r − m n + (1 − α)1 q, forα satisfying j/m ≤ α ≤ 1.

We can find the relation by matching the dimensions of the above inequatity. U Ln/p−j =U Ln/r−mα U Ln/q1−α= U L(n/r−m)α+(1−α)(n/q) This gives n p − j = α n r − m + (1 − α)n q.

1.5 Scaling

In the differential equations, we usually have some small parameters (or large parameters) after performing scaling. Roughly speaking, the high order derivative terms are important in a small spatial region or short time, while low order terms are important in large domain or long time. For instance, consider the following second order equation

y00+ ay0+ by = 0.

If we perform the scaling x∗ = x/. This means that in a small region of O() width, x is blowup to x/. The equation in terms of x∗is

1 2 d2y dx∗2 + 1 dy dx∗ + by = 0.

In this equation, the highest order term is the most important term. On the other hand, if we rescale x = λx0with λ >> 1, then the lowest order term is the most important term.

1.6 Homework

• pp. 17: 1, 2, 3 • pp. 19: 12, 15.

–

– The dimension of mass concentration is [C] = M L−3. • pp. 30: 3, 6, 7, 13

(18)

1.7. PROJECTS 15

1.7 Projects

Choose one of the following subject for your project. This should be turned in near the end of this semester.

1. Study Benard Flow and its dimensional analysis. Ref.

2. Study atmospherical flows and the corresponding dimensional analysis 3. Study Komogorov’s turbulence theory.

4. Study water wave theory and its dimensional analysis 5. Study phase field model for multiphase flows

6. Study Schr¨odinger equation and its semi-classical scaling.

7. Study the Ginzburg-Landau theory for superconductivity, perform dimensional analysis and derive a dimensionless form.

(19)

(20)

Chapter 2

Perturbation Methods

2.1 Perturbation methods for algebraic equations

Let us suppose our algebraic equations depend on a parameter . Suppose the root can be found for = 0. We look for roots for small . The procedure of regular perturbation are the follows:

• Express x = x0+ x1+ 2x2+ · · · ;

• Plug this expression into equation, make Taylor expansion of coefficients of the equation. • Equating the coefficients with like power of ;

• Solve x0, x1successively.

Example 1 Consider x(x − 1) = .

• Let x = x0+ x1+ · · · . Plug this into equation, we get

(x0+ x1+ · · · )(x0+ x1+ · · · − 1) = .

• Equating the coefficients of like powers:

0 : x0(x0− 1) = 0,

1 _: _2x

0x1− x1= 1,

.. . : ... This leads to two sets of solutions:

x(1)_a = 0 − ; x(2)_a = 1 + The true solution is

x = 1 ± √ 1 + 4 2 ≈ 1 2(1 ± (1 + 2)). which are consistent to the regular perturbation solutions.

(21)

18 CHAPTER 2. PERTURBATION METHODS Example 2 Let us consider the equation x2 = . You will see that

• The expansion x = x0+ x1+ · · · does not work. You should try x =

√ . • For < 0, the solution becomes imaginary.

So, we should try x = x0+ δ()x1+ δ()2x2+ · · · .

• Plug this ansatz, we get

x2₀+ 2x0δx1+ δ2x21+ 2δ2x0x2+ · · · = .

By comparing both sides, we see that we should have x0 = 0. Then this gives

δ2x2₁+ · · · =

To equate both sides, we need to choose δ2= . This give δ =√ and x1 = 1.

Example 3. Let us consider

x2+ ax + b = 0.

As → 0, there is only one root. Thus, the perturbation method can not recover the other root, which goes to ∞ as → 0. If a = 0, there is no root as → 0, the reduced equation is even inconsistent at all. Thus, the above perturbation method does not work for such case.

Nevertheless, we can try the following thing. We know the other solution goes to infinity as → 0, we try x =x∗. Plug this into equation, we obtain

x∗2 + a

x∗

+ b = 0.

We see that x∗satisfies an equation where the regular perturbation method can handle. We write x∗= x∗₀+ x∗₁+ · · · ,

Plug this into the above rescaled equation, we get

x∗₀2+ ax∗ = 0 2x∗₀x∗₁+ ax∗₁+ b = 0

These equations give the x∗₀ = −a and x∗₁= b/a, lead to the second solution of the orginal equation x(2)= −a

+ b

a+ · · · . You may think how to handle the case when a = 0.

(22)

2.2. REGULAR PERTURBATION METHOD FOR DIFFERENTIAL EQUATIONS 19 Homework

1. Find the asymptotic behaviors of the equations (a) x3+ x − 2 = 0;

(b) 2x4+ x3+ x − 1 = 0 (c) x5+ x3− 1 = 0;

2. There are infinite roots of tan x = x. Find their asymptotic fomula.

2.1.1 Justification of regular perturbation method for algebraic equations

Implicit Function Theorem

Theorem 2.3. Let F : Rn+1 → Rn _{be smooth. Suppose} _x

0 is a solution of F (x0, 0) = 0

and suppose ∂F/∂x(x0, 0) is non-singular. Then there is a smooth solution set x() satisfying

F (x(), ) = 0 for small .

The proof of this theorem is based on method of contraction map. We rewrite the above equation as a perturbtion equation: let us write x() = x0+ y(), then y() satisfies

∂F ∂x(x0, 0)y + ∂F ∂(x0, 0) + r(y) = 0. Here, r(y) = F (x0+ y, ) − F (x0, 0) − ∂F ∂x(x0, 0)y − ∂F ∂(x0, 0) = O(|y| 2₊2_).

Since J := ∂F/∂x(x0, 0) is non-singular, we take its inversion and get

y = T y := J−1 −∂F ∂(x0, 0) − r(y) .

We want to find a small number 0and another number η such that for any || ≤ 0, the mapping T

is a strict contraction map from |y| ≤ η to itself. Then by the fixed point theorem, we can obtain a fixed point y().

This process also tell us the construction of the perturbed solution. The method breaks down when the Jacobian J := ∂F_∂x(x0, 0) is singular, or when it has very small eigenvalue.

Extended study

1. A short note on regular perturbation method by Eric Vanden-Eijnden.

2.2 Regular perturbation method for differential equations

(23)

20 CHAPTER 2. PERTURBATION METHODS A falling object with resistivity The model reads

mdv

dt = −av + bv

2_{, v(0) = V} 0.

We introduce the dimensionless variables y = v/V0, τ = at/m, then the equation becomes

dy dτ = −y + y 2_, where = bV0 a << 1.

It means the resistivity (damping) is very large, as compared with V0and b. This equation has exact

solution

y(τ ) = e

−τ

1 + (e−τ_{− 1)}.

which has the following Taylor expansion in :

y = e−τ+ (e−τ− e−2τ) + 2(e−τ2e−2τ+ e−3τ) + · · · . The regular perturbation method introduces a Taylor expansion of y in terms of :

y(τ, ) = y0(τ ) + y1(τ ) + 2y2(τ ) + · · · .

We plug this ansatz into the equation, equating the coefficients of like powers of . We get y₀0 = −y₀,

y₁0 = −y₁+ y2₀ y₂0 = −y2+ 2y0y1,

.. . Equating the initial conditions, we get

y0(0) = 1, y1(0) = y2(0) = · · · = 0.

Solving these equations, we get

ya= e−τ+ (e−τ − e−2τ) + 2(e−τ2e−2τ + e−3τ) + · · · ,

We find this approach does work for this example.

The logistic model We consider the logistic model for population dynamics: dP

dt = P (1 − P K)

(24)

2.2. REGULAR PERTURBATION METHOD FOR DIFFERENTIAL EQUATIONS 21 Nonlinear oscillator Consider a nonlinear oscillator

md

2_y

dτ2 = −ky − ay

3_{, y(0) = A,}dy

dτ(0) = 0. This is so-called hard spring. We rescale it by

t = τ m/k, u =

y A. Then we get the Duffing equation:

d2u

dt2 + u + u 3 _{= 0}

u(0) = 1, u0(0) = 0. We perform regular perturbation method

u(t, ) = u0(t) + u1(t) + 2u2(t) + · · · ,

Plugging this into equation and the initial conditions, we get ¨

u0+ u0 = 0, u0(0) = 1, ˙u0(0) = 0,

¨

u1+ u1 = −u30, u1(0) = 0, ˙u1(0) = 0,

From the fisrt equationm, we obtain

u0(t) = cos t.

The second equation becomes ¨

u1+ u1 = − cos3t = −

1

4(3 cos t + cos 3t) Solving this equation with initial condition, we get

u1(t) =

1

32(cos 3t − cos t) − 3 8t sin t.

The term t sin t is a resonant term from cos t. Such a term is called a secular term. It will grow linearly and eventually to infinite as t → ∞. However, by energy method, one can show that the solution is bounded. What wrong is that the expansion is only good for finite time. The estimate |ya(t, ) − ye(t, )| = O(2) is only valid for t ∈ [0, T ] for a finite T .

Asymptotic Expansion First, we give some definitions of some notations. • The notation f () = o(g()) means that f /g → 0 as → 0.

• If f and g are also function of t in an interval I, the notation f (t, ) = o(g(t, )) as → 0, t ∈ I means that for every t ∈ I, we have f (t, )/g(t, ) → 0 as → 0.

(25)

22 CHAPTER 2. PERTURBATION METHODS • If the above limit is uniform for t ∈ I, we say that f (t, ) = o(g(t, )) as → 0 uniformly on

t ∈ I.

Definition 2.1. • A sequence of gauge functions {g_n(t, )} is an asymptotic sequence on t ∈ I as → 0 if

gn+1(t, ) = o(gn(t, )) for every t ∈ I.

• Given a function y(t, ) and an asymptotic sequence {g_n(t, ) on t ∈ I, the formal expansion P∞

n=0angn(t, ) is said to be an asymptotic expansion of y(t, ) as → 0 if

y(t, ) −

N

X

n=0

angn(t, ) = o(gN(t, )), as → 0,

for anyN . If the limits are uniform for t ∈ I, we say it is a uniform asymptotic expansion on I.

The expansion sequences are usually seperable such as • n_u

n(t)

• αn_u

n(x), where αnis a strictly increasing sequence.

• n_{ln u} n(x).

A rigorous proof for regular perturbation method has been done. Basically, the highest order term is elliptic operator on finite domain, the perturbation should be in the low orders and can be controlled by the elliptic operator. For instance, consider

4u = f (u, Du, )in Ω, u = 0 on ∂Ω.

We assume the equation is solvable for = 0. We also assume f (u, Du, ) can be controlled by 4u. This means that 4−1_{f (u, Du, ) is a compact smooth map from, say H}1 _{to H}1_{. Then we can}

apply implicit function theorem to get the solution for small . In general, the term Du is harder to control, it relies on Sobolev embedding theorem.

If the underlying operator is wave opertor, or Schr¨odinger, then it is even harder. There is no such compactness property. Nash-Moser technique is introduced.

2.3 The Poincar´e-Lindstedt Method

In the Duffin’s equation: ¨

u + u + u3= 0, u(0) = 1, ˙u(0) = 0,

the regular perturbation method leds to a secular term, which is incorrect for large time. One way to solve this problem is to introduce a change of time scale. Let

(26)

2.4. SINGULAR PERTURBATION METHODS 23 τ = ω()t, ω = 1 + ω1+ 2ω2+ · · · ,

Then

ω2u00+ u + u3 = 0, u(0) = 1, u0(0) = 0,

where0 = d/dτ . We plug the expansion above for u and ω into this equation, equating the coeffi-cients of likely powers of , we get

(1 + 2ω0ω1+ · · · )(u000+ u 00 1+ · · · ) + (u0+ u1+ · · · ) + (u30+ 3u20u1+ · · · ) = 0, u0(0) + u1(0) + · · · = 1, u00(0) + u01(0) + · · · = 0, u00₀+ u0= 0, u0(0) = 1, u00(0) = 0, · · · , u00₁ + u1 = −2ω1u000− u30, u1(0) = u01(0) = 0, This gives u0(τ ) = cos τ. u00₁+ u1 = 2ω1cos τ − cos3τ = (2ω1− 3 4) cos τ − 1 4cos 3τ. If we choose ω1= 3/8, then the secular term can be avoided. This leads

u1(τ ) =

1

32(cos 3τ − cos τ ). Thus, we get the expansion

u(τ ) = cos τ + 32(cos 3τ − cos τ ) + · · · , τ = t +3 8t + · · · . Homework. • pp. 101: 8(a), • pp. 102: 11, 13. 15.

2.4 Singular perturbation methods

2.4.1 Outer solutions, inner solutions and matched asymptotics

If the small parameter appears in the highest order term, then this term is unimportant in most of the region except in a small region where the high order derivatives are important. Let us see the following example:

(27)

24 CHAPTER 2. PERTURBATION METHODS u(0) = 0, u(1) = 1.

Here, we assume a < 0. Physically, is the viscosity, a the advection velocity. In our present situation, the advection direction is toward left. This equation can be solved easily. We integrate it once to get

−au + u_x = C1,

where C1 is a constant to be determined. Using method of separation of variable,

du C1+ au

= dx . Integrate this again, we get

1 aln(au + C1) = x + C2. This gives u = C3exp ax + C4.

Putting the boundary conditions, we get

C3+ C4 = 0, C3ea/+ C4= 1. These gives C3 = 1 ea/_{− 1}, C4 = − 1 ea/_{− 1}.

Hence, the exact solution is

ue(x) =

1 − eax/ 1 − ea/ .

Next, we shall use perturbation method to find approximate solution. It consists of three steps: (1) finding outer solution, (2) finding inner solution, (3) matching the outer and inner solution. Finding Outer solution Let us first solve this equation with = 0:

aux = 0.

This leads to u = constant. From boundary conditions, they are two possible solutions, u = 0 or u = 1. As we shall see later that the condition a < 0 leads to the flows move toward left and hence we should use the boundary condition u(1) = 1. This implies u(x) ≡ 1 for the unperturbed equation. This also means that u(x) ≡ 1 in the region where u is smooth (hence uxx is small. The

solution

uo(x) = 1

(28)

2.4. SINGULAR PERTURBATION METHODS 25 Finding inner solution However, u cannot be smooth through out the whole region (0, 1) because u(0) = 0. Therefore, we expect there is an abrupt change of u near x = 0. This region is called boundary layer. Let us suppose its thickness is δ(). We rescale ξ = x/δ. The rescaled equation becomes δ2 d2u dξ2 − a δ du dξ = 0.

If we want to have these two terms to be equally important, then we should take

δ2 =

1 δ. In this case, δ = . The resulting equation is

uξξ= auξ

with boundary condition u(0) = 0. Now, we get the same equation without . This gives uξ = au + C1.

The solution is

u(ξ) = Ceaξ− 1. where C is to be determined. The solution

ui(x) = C

eax/− 1

is called inner solution. This inner solution is an approximate in the region (0, ). To determine C, we need to match ui(x) and uo(x) for x in some overlapping zone. A natural overlapping zone is

when η = x/√ with η = O(1). In this case, let us fix η and we expect uo(

√ η) − ui( √ η) → 0 as → 0. This implies lim →0C eaη/ √ _{− 1}_{= 1.}

Since a < 0, we get from the above limit that C = −1. Thus, ui(x) = 1 − eax/.

Match inner and outer solutions We can define an approximate solution to be the sum of the inner solution and outer solution minus the overlapping value. The overlapping value is 1. Thus, we define

ua(x) = uo(x) + ui(x) − 1 = 1 − eax/.

We see that ua(0) = 0, ua(1) = 1 − ea/. In fact,

ue(x) − ua(x) = 1 − eax/ 1 1 − ea/ − 1 = 1 − eax/ e a/ 1 − ea/

(29)

26 CHAPTER 2. PERTURBATION METHODS Remark. Notice that the above approximate solution ua(x) does not satisfy the boundary

condi-tion at x = 1, but has a small error. It does satisfy the equacondi-tion and the boundary condicondi-tion at x = 0. Alternatively, we can choose different approximate solution. For instance, let ω(x) be a weighted function which is 1 for 0 ≤ x ≤ 1/3 and 0 for 2/3 ≤ x ≤ 1 and smoothly and monotonely connect 1 and 0 for 1/3 ≤ x ≤ 2/3. Using this weighted function, we define

ua(x) = ω x β ui(x) + 1 − ωx β uo(x).

Here, β is any number between 0 and 1. This approximate solution satisfies the boundary condition, but does not satisfy the equation, with a small residual.

Example This example is taken from J. Cole’s book, Singular Perturbation, pp. 21. Consider uxx+

√

xux− u = 0, x ∈ (0, 1),

u(0) = 0, u(1) = e2.

1. Finding outer solution: just like regular perturbation method, we assume u(x) = u0(x) + u1(x) + · · · .

Plugging this into equation and the boundary conditions, equating the coefficients of the like power terms of , we get _√

xu0,x = u0, u0(1) = e2,

√

xu1,x = u0− u0,xx, u1(1) = 0.

This leads to the outer expansion: uo(x) = e2 √ x 1 − 1 2x− 2 √ x + 3 2 + O 2 x5/2

The reason why we only use the boundary condition at x = 1 is because the advection velocity is negative (−√x), which means that the upwind direction is right. By the method of characteristics, the solution is determined by its upwind data, thus, from the right. Hence the boundary condition for the outer solution is from x = 1.

2. Finding inner solution: The boundary layer occurs near x = 0. We rescale x by introducing the layer variable

ξ = x δ(). The boundary layer expansion is

u = w0(ξ)ν0() + w1(ξ)ν1() + · · · ,

Plug this into the equation, we get

δ2 (ν0w0,ξξ+ ν1w1,ξξ+ · · · ) + √ δ δ √ ξ (ν0w0,ξ+ ν1w1,ξ + · · · ) −ν0w0− ν1w1− · · · = 0.

(30)

2.4. SINGULAR PERTURBATION METHODS 27 We choose δ() so that δ2 = √ δ δ . This gives δ() = 2/3. The dominant boundary equation becomes

w0,ξξ+ p ξw0,ξ = 0, w0(0) = 0. Its solution is w0(ξ) = C0 Z ξ 0 exp −2 3ζ 3/2 dζ. Thus, the inner solution has the form:

ui(x) = w0

x 2/3

+ · · · . 3. Matching: Let the matching scale is

xη =

x η η() is chosen so that, with fixed xη, as → 0,

• the outer variable x = ηxη → 0;

• the inner variable: ξ = x/δ = xηη/δ → ∞.

For instance, we can choose η = β with 0 < β < 2/3. The outer and inner solution should match in the overlapping zone where xηis fixed. As → 0,

• The outer solution u0(x) → 1, as x → 0.

• The inner solution

w0(ξ) → C0 Z ∞ 0 exp −2 3ζ 3/2 dζ, as ξ → ∞. To match these two limits, we should require

C0= Z ∞ 0 exp −2 3ζ 3/2 dζ −1 .

This gives the complete description of inner solution. The approximate solution is then de-fined to be

outer solution + inner solution − overlapped value That is, ua(x) = uo(x) + w0 x 2/3 − 1.

(31)

28 CHAPTER 2. PERTURBATION METHODS Remark. The next term is ν1() = 1/3. Check by yourself.

2.4.2 The boundary layers, initial layers and interior layers

Initial layer Let us consider a damped spring-mass system: m¨y + a ˙y + ky = 0,

with y(0) = 0 and m ˙y(0) = I. This means that we apply an impulse at the mass at t = 0. The dimensions of these variables are

[m] = M, [a] = M T−1, [k] = M T−2, [I] = M LT−1. Three possible time scales are

m a, r m k, a k,

corresponding to balancing inertia and damping, inertia and spring stiffness, damping and spring stiffness. Possible length scales are

I a, I √ mk, aI mk

We expect that the impulse will cause an abrupt change of the mass position in short time, then relax to its equilibrium. The first time period is called an initial phase, the second is called an relaxation phase. In the initial phase, the dominated terms should be the inertia term and the damping terms. In the relaxation phase, it should be a balancing between damping and spring stiffness. Thus, for the realxation phase, we introduce the time scale

¯ t = t

a/k. The equation becomes

ma2 k2 d2y d¯t2 + a2 k dy d¯t + ky = 0,

In the initial phase, the amptitude of the mass is related to I and the damping, and the mass too. However, from dimensional analysis, the mass M appears in both I and a. Thus, the amptitude should be related only to I/a. Thus, in the relaxation phase, we rescale the length by

¯ y = y

I/a. The rescaled equation becomes

¯y00+ ¯y0+ ¯y = 0, ¯

y(0) = 0, ¯y0(0) = 1. Here, the dimensionless parameter

= mk a2 << 1.

(32)

2.4. SINGULAR PERTURBATION METHODS 29 The outer solution is a solution of ¯y0+ ¯y = 0, this gives the outer solution

¯

yo(¯t) = Ce−¯t.

During the initial phase, we rescale τ = ¯t/ and Y = ¯y. Then d2Y

dτ2 +

dY

dτ + Y = 0. The conditions ¯y(0) = 0, ¯y0(0) = 1 gives the inner solution

¯ yi(¯t) = Y ( ¯ t ) = 1 − e −¯t/_.

Matching the outer solution and inner slution in an overlapping zone, (i.e. ¯tη = ¯t

√

is fixed), we should require

lim

→0y¯0(¯t) = lim→0y¯i(¯t).

This leads to C = 1 in the outer solution. Thus, the final approximate solution is ¯

ya(¯t) = y¯o(¯t) + ¯yi(¯t) − lim →0y¯i(¯t)

= e−¯t− e−¯t/. In terms of the original variables, it reads

ya(t) =

I a

e−kt/a− e−at/m.

Enzyme Kinetics Consider the following chemical reaction: A + B C → P + B

Here, A (a substrate) and B (an enzyme) combine to form a molecule C. C breaks into a product P and an orginal enzyme B. Let a, b, etc represent their concentration. The equations for the kinetics are ˙a = −k₁ab + k2c ˙b = −k1ab + k2c + k3c ˙c = k1ab − k2c − k3c ˙ p = k3c.

The initial concentrations are

(33)

30 CHAPTER 2. PERTURBATION METHODS It is easy to see that b + c = ˆb. The first three equations do not involve p. Thus, we have the reduced equations

˙a = −k1a(ˆb − c) + k2c,

˙c = k1a(ˆb − c) − (k2+ k3)c.

The dimensions of each quantities are

[a] = C, [c] = C, [k1] = T−1C−1, [k2] = [k3] = T−1.

Here, C is the dimension of concentration. We can rescale a and c by ¯a = a/ˆa and ¯c = c/ˆb. We also rescale time by ¯t = t/T . The dimensionless equation becomes

ˆ a T d¯a d¯t = −k1aˆˆb¯a(1 − ¯c) + k2ˆb¯c, ˆ b T d¯c d¯t = k1ˆaˆb¯a(1 − ¯c) − (k2+ k3)ˆb¯c,

To determine the time scale T , we notice that our interest is to study how A is converted to P through the enzyme B. Thus, the decreasing time scale should be the time scale we should concern. Thus, we should balance_Taˆ and k1ˆaˆb in the equation for ¯a. Thus, T is taken to be

T = 1 k1ˆb

. With this time scale, we obtain the dimensionless equations:

d¯a d¯t = −¯a + (¯a + λ)¯c, d¯c d¯t = ¯a − (¯a + µ)¯c. Here, = ˆb ˆ a << 1, λ = k2 ˆ ak1 , µ = k2+ k3 ˆ ak1 . The initial conditions are

¯

a(0) = 1, ¯c(0) = 0.

The outer solution is obtained by setting = 0. From this, we obtain ¯a − (¯a + µ)¯c = 0, which gives ¯

c = a¯ ¯ a + µ, and the first equation becomes

d¯a d¯t = −

µ − λ 1 +µ_¯_a

(34)

2.4. SINGULAR PERTURBATION METHODS 31 By separation of variable and integrating it, we get

¯

a + µ ln ¯a = −(µ + λ)¯t + K.

Here, the constant K will be determined from matiching with the inner solution. These two solutions a and c are our outer solutions. We denote them by aoand corespectively.

In the initial layer, we use the non-dimensional variables τ = t

, A = a, C = c. The resulting equations

dA

dτ = (−A + (A + λ)C) dC

dτ = A − (A + µ)C

Setting = 0, we obtain dA/dτ = 0. From the initial condition ¯a = 1, we should take A = 1. Plug this into the second equation, we get

dC

dτ = 1 − (1 + µ)C. With the initial condition c(0) = 0, we obtain

C = 1 − e

−(µ+1)τ

µ + 1 .

For matching, we should have the outer solution ao(0) = A(∞) and co(0) = C(∞). After some

calculation, we get K = 1 from matching condition. Thus, the final approximate solution is aa(t) = ao(¯t)

ca(t) = co(¯t) + C(τ ) −

1 µ + 1.

Boundary layers and internal layers We have seen that the Sturm-Liouville system uxx− a(x)ux− b(x)u = 0 x ∈ (0, 1),

u(0) = u0, u(1) = u1,

may have boundary layer. We will show that

• if a(x) > 0 for x ∈ (0, 1), then the boundary layer appears at x = 1; • if a(x) < 0 for x ∈ (0, 1), then the boundary layer appears at x = 0;

• if a(x) changes sign once in (0, 1) with a(0) > 0 and a(1) < 0, then an internal layer is formed.

(35)

32 CHAPTER 2. PERTURBATION METHODS Shock wave Consider the Burgers’ equation

uux= uxx, x ∈ (−1, 1)

with u(−1) = 1 and u(1) = −1. We can see that the outer solution is uo(x) =

1 −1 ≤ x < x₀ −1 x0 < x ≤ 1

The constant x0will be determined later.

For inner solution, we rescale it by ¯x = (x − x0)/. The equation becomes

u2 2 ¯ x = ux¯¯x.

We integrate it once to get

u2

2 − ux¯ = C.

Applying matching condition (u(±∞) = ∓1), we obtain C = 1/2. Using separation of variable, we get 2du u2_{− 1} = d¯x. Integrating ln u + 1 u − 1 = ¯x − ¯x0

Here, ¯x0is a constant. We can absorb ¯x0into x0. Thus, we take ¯x0= 0. Then

u + 1 u − 1 = e ¯ x_. This yields u(¯x) = 1 + e ¯ x 1 − ex¯ = − tanh(¯x).

The parameter x0 is not unique. Any such interior layer solution with the outer solution forms an

approximate solution unless we impose an extra condition. Usually, we impose an excess mass based on conservation of mass. This means that

Z 1

−1

u(x) − uo(x) dx = m

where m is called the excess mass. With this, ¯x0is determined by

Z x0 −1 − tanh x − x0 − 1, dx + Z 1 x0 − tanh x − x0 + 1, dx = m By taking → 0, we can obtain an approximation for x0.

x0 =

m

(36)

2.5. WKB METHOD 33 Homework The Burgers equation:

ut+ uux = uxx

is a prototype in continuum mechanics, turbulence theory, etc. Here we consider the steady Burgers’ equation

uux= uxx, x ∈ (0, 1),

with the boundary condition

u(0) = u0, u(1) = u1

Find its asymptotic solutions. You need to classify different cases according to the signs of u0and

u1.

2.5 WKB method

The WKB method is a perturbation method for solving problems of the followng form −2u00+ q(x)u = 0

When q < 0, we expect exponential decay solution. When q > 0, we expect oscillatory solution. In both cases, we look solution of the form ew(x).

Nonoscillatory Case Consider

2u00− k(x)2u = 0, x ∈ (0, ∞) Let us try regular perturbation for large x. We write

u = u0+ u1+ 2u2+ O(3),

Plug into the equation, we get −k2_u

0− k2u1+ 2(u000− k2u2) + · · · = 0,

This leads to u0 = u1 = u2 = · · · = 0. We get no information from regular perturbation for large

x. If we observe the equation more carefully, suppose k > 0, then we expect the solution decays for large x. Thus, let us figure how it decays by trying the ansatz u = ew(x). Then

u0 = w0ew, u00= (w00+ w02)ew, and the equation becomes

2(w00+ w02) − k2ew = 0. Thus, as long as w 6= −∞, we get

(37)

34 CHAPTER 2. PERTURBATION METHODS Let us introduce v = w0We have

v0+ v2− k(x)2= 0. Apply the regular perturbation method

v = v0+ v1+ · · · We get v0(x) = ±k(x), v1(x) = − k0 2k. This gives w0= ±k(x) − k 0 2k + O( 2₎ Or w(x) = 1 ± Z x a k(ξ) dξ − 2ln k(x) + O( 2₎ . Thus, we get an expansion

u(x) = 1 pk(x)exp ±1 Z x a k(ξ) dξ + O() = 1 pk(x)exp ±1 Z x a k(ξ) dξ (1 + O()) .

If we require u(∞) to be bounded, then we can only accept the exponential term. Suppose k > 0, for instance, then

u(x) = C 1 pk(x)exp −1 Z x a k(ξ) dξ (1 + O()) is admissible.

For x ∼ 0, we can rescale x by x0 = x/. Expand u(x0) in Taylor series near x0 = 0. We get −2u₂− 6u₃x0+ (k0+ k1x0+ ...)(u0+ u1x + ...) = 0.

This leads to

k0u0− 2u2 = 0

−6u₃+ k0u1+ k1u0 = 0,

We can determine u0 from the boundary condition. We can determine u1 from the matching with

(38)

2.5. WKB METHOD 35 Oscillatory cases For the Schr¨odinger equation

2u00+ k(x)2u = 0, where k > 0, we can approximate u by eiw(x). Then

u0 = iw0eiw, u00= (iw00− w02)eiw, and equation becomes

2(iw00− w02) + k2 = 0. If we introduce v = w0, then

iv0− v2+ k2= 0. The ansatz for the regular perturbation v = v0+ v1+ · · · gives

v0 = ±k, v1= iv₀0 2v0 = i(ln√v0)0. w(x) = ±1 Z x a k(ξ) dξ + i lnp±k(x) + O(). u(x) = exp ±i Z x a k(ξ) dξ exp − lnp±k(x)(1 + O()) . Thus, u(x) = c1 pk(x)exp i Z x a k(ξ) dξ + c2 pk(x)exp −i Z x a k(ξ) dξ .

2.5.1 Method of geometric optics

In optics, the governing equation is

utt= c(x)24 u.

For waves with a fixed frequency ω, the solution has the form: u(x, t) = e−iωtu(x) and u(x) satisfies the Helmholtz equation

4u + ω

2

c(x)2u = 0.

We are interested in the high frequency approximation of its solution. Let us rewrite ω = 1/. The ansatz is

u(x) = A(x, )eiφ(x)/ = (A0(x) + A1(x) + · · · )eiφ(x)/

where φ is a phase function, A the amptitude. We plug u into the Helmholtz equation to get uxi = Axie

iφ/₊iA

φxie

(39)

36 CHAPTER 2. PERTURBATION METHODS 4u = 4A +2i ∇A · ∇φ + iA 4 φ − A 2|∇φ| 2 eiφ/ −A 2|∇φ| 2₊ i (A 4 φ + 2∇A · φ) + 4A + A 2_c(x)2 = 0.

Expanding A in A0+ A1+ · · · and equating the coefficients of , we get

|∇φ|2 ₌ 1

c(x)2,

A04 φ + 2∇A0· ∇φ = 0.

The first equation is called the eikonal equation. The second equation is called the transport equa-tion. In the second equation, if we rename A0 =

√

ρ0, then the transport equation becomes

√ ρ04 φ + 2∇ √ ρ0· ∇φ = √ ρ04 φ + 1 √ ρ0 ∇ρ₀· ∇φ = 0. Thus, this is equivalent to

∇ · (ρ0∇φ) = 0.

WKB method for Schr¨odinger equation Consider the Schr¨odinger equation i~∂tψ = −~

2

2m∇

2_{ψ + V (x)ψ.}

We look solution of the form: ψ = ew/~. We have ∂tψ = 1 ~∂twe w/~_, ∇ψ = 1 ~ ∇wew/~, ∇2ψ = 1 ~ ∇2w + 1 ~2 X k (∂xkw) 2 ! ew/~. Thus, we get i∂tw = − 1 2m ~∇ 2_{w +}X k (∂xkw) 2 ! + V. We write w = R + iS, real part and imaginary part. Then

i∂t(R + iS) = −

1 2m ~∇

2_{(R + iS) + |∇R|}2_{− |∇S|}2_{+ 2i∇R · ∇S + V.}

Equating the real part and imaginary part, we get ∂tR = −

1 2m ~∇

(40)

2.5. WKB METHOD 37

−∂_tS = − 1 2m ~∇

2_{R + |∇R|}2_{− |∇S|}2_{+ V.}

1 _{Since [S] = ET = M L}2_T−1_{, [∇S] = M LT}−1_{, we can define p = ∇S, v = (∇S)/m.}

Multiplying the first equation by e2R/~, we get e2R/~∂tR = − 1 2m ~e2R/~∇2S + 2e2R/~∇R · ∇S = − ~ 2m∇ · e2R/~∇S

Next, we define ρ = e2R/~, use v = ∇S/m, then the first equation can be written as ∂tρ + ∇ · (ρv) = 0.

This is the continuity equation. For the second equation

St+ 1 2m|∇S| 2_{+ V =} 1 2m(|∇R| 2 + ~∇2R) We express the RHS in terms of ρ. We get

St+ 1 2m|∇S| 2_{+ V =} ~2 4m ∇2_ρ ρ . We take gradient of this equation to get

mvt+ m∇ |v|2 2 + ∇V = ~ 2 4m∇ ∇2_ρ ρ vt+ v · ∇v + 1 m∇V = ~2 4m2∇ ∇2_ρ ρ

Multiplying ρ to this equation, we get ρ(vt+ v · ∇v) + ρ m∇V = ρ~2 4m2∇ ∇2_ρ ρ

The left-hand side is the pressureless momentum equation for the gas dynamics. The right-hand side is a dispersion term. It regularizes the solution.

Thus, solving the outer solution based on the WKB approach for the Schr¨odinger equation is equivalent to solve a pressureless Euler equation. Although it is also not an easy job, the WKB approach provides us an insight of the macroscopic structure of the Schr¨odinger equationm. It also make a link between qnautum mechanics and classical mechanics.

1

The dimension of w and ~ are action, which energy times time ET . The dimension of the first equation is [R]M T−1= [~][S]L−2= [R][S]L−2, Thus, this is consistent to E = T−2L2M.

(41)

38 CHAPTER 2. PERTURBATION METHODS Reference

1. Bender-Orszag’s book has a complete description of WKB method for Sturm-Liouville equa-tion (in one dimension). In particular, it contains the descriptio of turning point.

2. Griffith’s book on Quantum mechanics has a good description on application of WKB for Schr¨odinger in one dimension.

Homework

• pp. 141: 2, 3, 5, 6,8.

2.6 Asymptotic expansion of integrals

2.6.1 Laplace method for approximation of integrals

In applied mathematics, many solution or approximate solution has integral representation. For instance, using WKB method, we can have approximate Green’s function formula (See Bender-Orszag). Using this formulation, we can have an integral representation of the solution. In particular, the following Laplace integral appears quite often in one dimensional problem:

I(λ) = Z b

a

f (t)e−λg(t)dt, λ >> 1. (2.1) And we are interested in its behavior for large λ. A generalization is when t is in the complex plane:

I(λ) = Z

C

f (z)eλg(z)dz, λ >> 1,

where C is a contour on complex plane. This leads to method of stationary phase and method of saddle point.

Laplace integrals We will study the asymptotic behavior of the integral of the form (2.1). The first observation is that the integrand e−λg(t) is important at those minima of g. So, let us first consider the case where g is monotonic in [a, b]. We may assume g is monotonic increasing in [a, b] without loss of generality. In this case, we can normalize the integral by the change of variable: s = g(t) − g(a). Then the integral becomes

I(λ) = e−λg(a) Z g(b)−g(a) 0 f (t(s)) g0(t(s))e −λs ds.

This is a standard form of the Laplace integral. The important part of the integral is near s 0. So, we decompose the integral into two regions: [0, T ] and [T, g(b) − g(a)]. In the latter region, the term e−λswill be exponential small (i.e. e−λT. For the former region, we make Taylor expansion of f (t(s))/g0(t(s)) near s 0, then integrate term by term to get expansion. We illustrate this approach by the following examples.

(42)

2.6. ASYMPTOTIC EXPANSION OF INTEGRALS 39 1. Example 1. Consider the integral

I(λ) = Z ∞ 0 sin t t e −λt dt, λ >> 1.

We decomposeR₀∞=R₀T +R_T∞. For the latter integral, we have Z ∞ T sin t t e −λt dt ≤ Z ∞ T e−λtdt = 1 λe −λT

For any fixed T , this term is always O(λ−N) for any N . We call it exponential small term (EST).

For the first integral, we make Taylor expansion of sin t near t = 0, we get

I(λ) = Z T 0 sin t t e −λt_{dt + EST} = Z T 0 N X n=0 (−1)n t 2n (2n + 1)! + O(t 2N +2₎ ! e−λtdt + EST = Z λT 0 N X n=0 (−1)n u 2n (2n + 1)!λ2n+1e −u_{du + O(λ}−2N −3₎ = Z ∞ 0 N X n=0 (−1)n u 2n (2n + 1)!λ2n+1e −u du + O(λ−2N −3) + EST = N X n=0 (−1)n Γ(2n + 1) (2n + 1)!λ2n+1 + O(λ −2N −3₎ = N X n=0 (−1)n (2n)! (2n + 1)!λ2n+1 + O(λ −2N −3₎ = N X n=0 (−1)n (2n + 1)λ2n+1 + O(λ −2N −3 )

Here, the function Γ(x) is called the Gamma function and is defined by Γ(x) :=

Z ∞

0

tx−1e−tdt.

This function is well defined for x > 0. By using integration by part, we get Γ(x + 1) = xΓ(x).

(43)

40 CHAPTER 2. PERTURBATION METHODS 2. Asymptotic expansion for error function. The error function is defined as

erfc(λ) = √2 π

Z ∞

λ

e−t2dt. We make a change of variable τ = (t − λ)/2. Then

erfc(λ) = √1 πe

−λ2Z ∞

0

e−τ2/4e−λτdτ.

We find the important part is localated near τ = 0. We can decompose the integral region into [0, T ] and [T, ∞). The latter part is a EST. In the formal integral, we make Taylor expansion for e−τ2/2for τ ∈ [0, T ], then integrate the integral, we get

erfc(λ) ∼ √2 πe −λ2 ∞ X n=0 (−1)n (2n)! n!(2λ)2n+1

Method of steepest descent for integrals In the case that g(t) has interior local minima, the important parts of the integral (2.1) are near these local minima. Let us suppose there is only one inimum, say t0. If there are more than one (but still finite many), we just add the asymptotic

formulation for each one together. We break the integral into (t0− δ, , t0+ δ) and its exterior part.

In the exterior region: |t − t0| ≥≥ δ and t is bounded,

g(t) − g(t0) ≥ M,

and the integral

Z

(a,b)∩|t−t0|≥δ

f (t)e−λt≤ Ce−λM.

For unbounded region (a, b), we can require that that the asymptotic behavior of g is some power growth, say g(t) − g(t0) ≥ M |t − t0|. Then we still have the above estimate.

In the interior region: |t − t0| < δ, let us expand g by

g(t) = g(t0) +

1 2g

00_(t

0)(t − t0)2+ O(t − t0)3, for |t − t0| ≤ δ.

We can replace g − g(t0) by the quadratic term in the Laplace integral:

Z |t−t0|<δ f (t)e−λg(t)= e−λ(g(t0)+O(δ3)) Z |t−t0|<δ f (t)e−λ(t−t0)2/2_dt. The term e−λO(δ3) = 1 + O(λδ3). The choice of δ should depend on λ. We shall determine it later.

For the integral

Z t0+δ

t0−δ

(44)

2.6. ASYMPTOTIC EXPANSION OF INTEGRALS 41 we make a change of variablepλg00_{(0)t − t}

0= s. Then the above becomes

Z δ √ λg00_(t 0) −δ√λg00_(t 0) f t0+ s pλg00_(t 0) ! e−s2/2d s pλg00_(t 0) ! .

We need the domain of integral goes to the whole line. This require δ√λ → ∞ as λ → ∞. On the other hand, we should also require λδ3 → 0 as λ → ∞. So, we can choose δ = λ−α _with

1/3 < α < 1/2. We may choose α = 1/2 − /3 for any small > 0. Then O(λδ3) = O(λ−1/2+).

With this, we get

I(λ) ∼ f (t0)e−λg(t0) 2π λg00_(t 0) 1/2 (1 + O(λ−1/2+)). Here, we have used

Z ∞

−∞

e−s2/2ds =√2π.

Remark. The error O(λ−1/2+) can be improved. Indeed, we can expand g as g(t) = 3 X i=0 g(i)(t0) i! (t − t − t0) i_{+ O(δ}4_).

We claim that the third order term can be removed by a qradratic translation. For simplicity, let us assume t0 = 0. We have g(t) = g(0) +g 00₍₀₎ 2 t 2₊g000(0) 3! t 3_{+ O(t}4₎ = g(0) +g 00₍₀₎ 2 t + at 22 + O(t4) where a = g 000₍₀₎ 3!g00₍₀₎.

Using the change-of-variable (t + at2) ← t, we can remove the third order term. Thus we get I(λ) = Z δ √ λg00_(t 0) −δ√λg00_(t 0) f (t0+ s pλg00_(t 0) )e−s2/2d s pλg00_(t 0) 1 + O(λδ4) . We can choose δ = M/√λ for large M . Then we can get the asymptotic formula

I(λ) ∼ f (t0)e−λg(t0) 2π λg00_(t 0) 1/2 (1 + O(λ−1)).

(45)

42 CHAPTER 2. PERTURBATION METHODS Stirning formula The Stirnling formula is

n! ∼√2πnn+1/2e−n(1 + O(1 n)).

This formula is very useful to connect discrete phenomena to continuous phenomena. Since Γ(n + 1) = n!, we shall study the asymptotic expansion of the Gamma function. The Gamma function is defined to be Γ(s) = Z ∞ 0 xs−1e−xdx = Z ∞ 0 e−sF (x,s)dx where F (x, s) = − ln x +ln x s + x s.

The minimum of F w.r.t. x is x = s − 1. We rescale x by setting t = x/(s − 1). Then F (t(s − 1), s) = − ln(t(s − 1)) + 1 sln(t(s − 1)) + t(s − 1) s = − ln(s − 1) + ln(s − 1) s + (− ln t + ln t s + t(s − 1) s ) Therefore, Γ(s) = (s − 1)sJ (s), where J (s) = Z ∞ 0 e−(s−1)g(t)dt, g(t) = t − ln t.

The minimum of g(t) is g0(t) = 1 − 1/t = 0. This gives the extremal point t0 = 1. At t = 1,

g00(1) = 1 > 0. Thus, t = 1 is a mimimum. From the Laplace asymptotic formula, we get J (s) ∼ e−(s−1) 2π s − 1 1/2 Thus, ln Γ(s) = s − 1 2 ln(s − 1) − (s − 1) + 1 2ln(2π) + O( 1 s − 1). Homework. • pp. 149: 3, 6, 7(a), 8, 11, 12.

(46)

Chapter 3

Calculus of Variation

3.1 Examples of Functionals

Many physical problems can be formulated as variational problems.

Example 1. Minimum arc length problem Consider the set

A = {y : [a, b] → R in C1, y(a) = ya, y(b) = yb}

For y ∈ A, we define its arc length by

J (y) = Z b

a

p

1 + y0_(x)2_dx

The function J maps a “path” y ∈ A to a real number. That is, J : A → R. We call J is a functional. The problem is to find a minimum of J in the set A. Certainly, the answer of this is the straight line connecting (a, ya) and (b, yb).

Example 2. Minimum area problem Consider the minimum area problem defined on a domain Ω ⊂ R2. Let

A = {u : Ω → R in C1, with u = u0, on ∂Ω}

where u0is a precribed function defined on ∂Ω. For a function u ∈ A, the area of its graph on Ω is

J (u) = Z

Ω

p

1 + |∇u(x)|2_dx.

The minimum of J in A is called the minimum surface with prescribed boundary value u0. This

problem is much nontrivial.

(47)

44 CHAPTER 3. CALCULUS OF VARIATION Example 3. The Brachistochrone The Brachistochrone problem is to find a curve on which a ball sliding down under gravitation to a point with depth h takes least time. The word “brachistochrone” means the “the shortest time delay” in Greek. It was one of the oldest problem in Calculus of Variation. Its solution is a section of a cycloid. This was founded by Leibnitz, L’Hospital, Newton and two Bernoullis. Suppose the curve starts from A = (0, 0) and ends at B = (a, −h). Let s be the arc length of the curve. The equation of motion is

m¨s = −mgy0(s). This gives the conservation of energy

1 2m ˙s

2_{+ mgy(s) = E.}

At point A, we take s = 0, ˙s = 0 and y(0) = 0. With this normalization, E = 0. Thus, we have the speed

v = ˙s =p−2gy.

Notice that y ≤ 0 under our consideration. The traveling time from A to B is given by

T (y) = Z T 0 dt = Z S 0 ds ˙s = Z s 0 ds √ −2gy.

Here, S is the arc length of the curve {(x, y(x))|x ∈ (0, a)}. We can also express the functional using parameter x: T (y) = Z S 0 ds v = Z a 0 p1 + y0_(x)2 p−2gy(x) dx The admissible class is

A = {y : (0, a) → R ∈ C1|y(0) = 0, y(a) = −h}. In the above examples, we encounter a functional defined as

J (y) = Z b

a

L(x, y(x), y0(x)) dx

for functions y in an addmissible class A. The class of functions should be continuously differ-entiable so that L(x, y(x), y0(x)) is well-defined. In addition, some boundary conditions of y at a and b are imposed. The integrand is called the Lagrangian. The functional J has the names “cost functional”, “action”, ... in different applications.

(48)

3.2. BANACH SPACES 45

3.2 Banach Spaces

Normed linear space The above functionals are defined in a function space. In abstract, the function space is a vector space V over a scalar field. The scalar field we consider is either the real field R or the complex field C. In a vector space, we can associat each vector y a norm kyk, which measures the length of the vector y. It is required to satisfy the following properties:

1. kyk ≥ 0, and kyk = 0 if and only if y = 0, 2. kαyk = |α|kyk,

3. Triangle inequality: ky1+ y2k ≤ ky1k + ky2k.

A vector space V associated with a norm k · k is called a normed linear space. Examples

1. Rnwith kxk =px2₁+ · · · + x2

nis a normed linear space.

2. Rnwith kxkp := (|x1|p+ · · · + |xn|p)1/pfor p ≥ 1 is a normed linear space.

3. C[a, b] = {y : [a, b] → R is continuous} with the norm kyk∞:= max

x∈[a,b]

|y(x)| is a normed linear space.

4. C1[a, b] := {y : [a, b] → R is continuously differentiable} with the norm kyk := max

x∈[a,b] |y(x)| + |y 0

(x)| is a normed linear space.

5. The space `p is defined to be

`p = {(xn)∞n=1|

X

n

|x_n|p _{< ∞}}

associated with the norm

kxkp := ∞ X n=1 |xn|p !1/p

The space `p with p ≥ 1 is a normed linear space. 6. The `∞space is

`∞= {(xn)∞n=1| sup n

|x_n| < ∞} associated with the norm

kxk∞:= sup n

(49)

46 CHAPTER 3. CALCULUS OF VARIATION 7. The space Lp(a, b) defined to be

Lp(a, b) = {y : (a, b) → R| Z b

a

|y(x)|pdx < ∞} associated with the norm

kykp :=

Z b a

|y(x)|pdx 1/p is a normed linear space for p ≥ 1.

8. The space L∞(a, b) is defined as

Lp(a, b) = {y : (a, b) → R| sup

x∈(a,b)

|y(x)| < ∞} associated with the norm

kyk∞:= sup x∈(a,b)

|y(x)| < ∞. 9. The Sobolev space Hm(a, b) is defined as

Hm(a, b) := y : (a, b) → R Z b a |y(x)|2+ · · · + |y(m)|2dx < ∞

with the norm

kyk_Hm := Z b a |y(x)|2_{+ · · · + |y}(m)_|2_dx 1/2 Homework 1. Show that kxkp → kxk∞as p → ∞.

Banach spaces In a normed linear space V , we can define the concept of limits.

• Limit: A sequence {yn} in V is said to have a limit y ∈ V if kyn− yk → 0. We denote it by

limn→∞yn= y.

• Cauchy sequence: A sequence {yn} in is called Cauchy sequence if all but finite of them are

closed to each other. In mathematical languish, for any > 0, there exists an N such that for all n, m > N , we have kyn− ymk < .

• Completeness: A normed linear space V is said to be complete if all its Cauchy sequences have limits in V . Such spaces are also called Banach spaces.

• Theorem [Completion of normed linear space] For any normed linear space V , there exists a minimal complete normed linear space ˜V containing V . One way to construct such a space is to define

˜

V = {(yn)|(yn) is a Cauchy sequence in V }/ ∼ .

Here ∼ is the notation of equivalence. Two Cauchy sequences (yn) and (zn) are called

(50)

3.2. BANACH SPACES 47 Examples

1. The space C[a, b] under the k · k∞is a Banach space.

2. The space C[−1, 1] under k · k1 is not a Banach space. Because the sequence tanh(x/)

with = 1/n is a Cauchy sequence, but its limit is a discontinuous function which is not in C[−1, 1].

3. The spaces Lp(a, b) under k · kp, 1 ≤ p ≤ ∞ are Banach spaces. In fact, the Lp(a, b),

1 ≤ p < ∞ is the completion of C[a, b] under the norm k · kp.

4. The Sobolev space Hm(a, b) the completion of Cm_{[a, b] under the norm:}

kykHm :=

Z b

a

|y(x)|2+ · · · + |y(m)|2dx 1/2

Homework: Find the function classes (as shown in the above examples) of the following functions they belong to.

1. The Heviside function (i.e. the step function).

2. Find the condition on α such that the function |x|α∈ Lp_{(−1, 1)}

Approximation and Subspaces In applications, we usually approximate a function by nice smooth functions. For instance, u ∈ C[a, b] can be approximated by C∞[a, b] by the norm k · k∞. That is,

for any > 0, we can find a function u ∈ C∞[a, b] such that ku − uk < . We can choose, for

instance, u= ρ∗ u, where ρ(x) := ρ(x/)/ and ρ ∈ C∞(R), has compact supported in (−1, 1),

ρ(x) > 0 for |x| < 1/2 andR ρ dx = 1.

The space C∞ ⊂ C[a, b] as a subspace. Moreover, from the discussion above, we say that C∞[a, b] is dense in C[a, b]. Some commonly used examples are:

1. Hm(a, b) is dense in L2(a, b) for 0 ≤ m < ∞,

2. Lp(a, b) is dense in Lq_{(a, b) for 1 ≤ q ≤ p ≤ ∞. Here, we require the interval (a, b) being}

bounded.

Homework. Prove the assertion: 1. u∈ C∞.

(51)

48 CHAPTER 3. CALCULUS OF VARIATION Boundary conditions and Subspaces In variational problems, we usually impose some boundary condition such as Dirichlet, Neumann or Robin conditions. They come naturally from the process of integration-by-part. We will explain this latter. We would like to point out here is that the set which satisfies the boundary condition usually form a subspace, provided the problem is linear. Here are some examples:

1. C0[a, b] = {u ∈ C[a, b], u(a) = u(b) = 0} is a closed subspace in C[a, b].

2. H₀1[a, b] = {u ∈ H1(a, b), u(a) = u(b) = 0} is a closed subspace in H1(a, b). 3. C₀m[a, b] = {u ∈ Cm[a, b], u(a) = u(b) = 0} is a subspace in Cm[a, b].

In general, the boundary condition gives a constraint. The set A of those vectors satisfying the constraint forms a manifold (or an admissible class) in the Banach space.

3.3 Linear operators

Linear functionals Let V be a normed linear space over R. Let ` be a linear function from V to R. Such a function is called a linear functional. A linear functional ` is said to be bounded if there exist a constant C such that |`(x)| ≤ Ckxk for all x ∈ V .

Proposition 1. ` is bounded if and only if it is continuous. Examples

1. Let x0 ∈ [a, b]. Consider the mapping defined by `(f ) = f (x0) for any f ∈ C[a, b]. Then `

is a bounded linear functional on C[a, b]. 2. Let > 0 be small number. Define

`(f ) = 1 2

Z

−

f (x0+ y) dy

Then ` is a bounded linear functional on L1(R).

3. Let ρ ∈ C∞with the properties: ρ > 0 in (−1/2, 1/2), supp ρ ⊂ (−1, 1), andR ρ = 1. Let be a small positive number, and define ρ(x) := 1ρ(x/). Define

`(f ) = ρ∗ f (x0) :=

Z

f (x0− y)ρ(y) dy.

Then ` is a bounded linear functional on L1(R).

4. The mapping: `(f ) := Df (x0) is a linear functional on C1[a, b] but not bounded. Here,

Methods in Applied Mathematics