Download Solution to Diffusion Equation , Lecture Notes - Advanced Calculus and more Study notes Calculus in PDF only on Docsity!
Solutions to the Diffusion Equation
Adrian Down
May 02, 2006
1 Review
Last time, we began to investigate the diffusion PDE,
∂tc = D∂xxc
Note. • ∂t represents (^) ∂t∂ , and ∂xx = (^) ∂x∂ 2.
- D is a constant that depends on the material in which the diffusion is taking place. We saw last time, by dimensional analysis, that D has units of, [D] = length
2 time. We considered the particular example of an initially point-like concentra- tion diffusing in water. To see the physical effect of diffusion, consider the case of sodium ions in water. In this case, the constant D takes the value D ≈ 1. 1 × 10 −5 cm
2 s. Using dimensional analysis, we saw that the size of the concentration grows like r(t) =
Dt. Taking t = 3 hours as an example,
r ≈
10 −^5
cm^2 s
· 104 s ≈ .3 cm
Hence diffusion is often a very slow process in practice.
2 Solving the diffusion equation
2.1 Motivation
It was mentioned at the end of the last lecture that the solutions to the diffusion equation take the form,
c(x, t) =
πDt
c(x′, 0)e−^
(x−x′)^2 4 Dt (^) dx′
The general strategy is to consider the evolution of the amount of sub- stance within a particular interval. Consider an initial amount of substance in the interval (x′, x′^ + dx′) at time t = 0. At time t > 0, this amount of substance will be distributed over a larger range in x′. In particular, the distribution of density takes the form,
density ρ ∝ e−^
(x−x′)^2 4 Dt
Because the material is spreading over time, the magnitude of the density is decreasing, which leads to a normalization factor,
ρ =
e−^
(x−x′)^2 4 Dt 2
πDt
The total amount of the substance is given by the product of the density and the initial amount of material in the interval (x′, x′^ + dx′),
dc(x′, t) dx′^
e−^
(x−x′)^2 4 Dt 2
πDt
c(x′, 0)
With this formulation, we see that the solution proposed above is the sum of the time-evolved distributions resulting from each small interval of the initial distribution.
2.2 Wave solutions
Previously, we solved recursion relations such as,
y′′^ + y′^ + y = 0
3.2 Finite difference approximations
We can approximate the derivatives of c using finite difference approxima- tions. We obtained these approximations when we studied Taylor series,
(∂xxc)(x, t) =
c(x + ∆x, t) − 2 c(x, t) + c(x − ∆x, t) (∆x)^2
+ O
(∆x)^2
We also found an approximation to the first derivative,
(∂tc) (x, t) =
c(x, t + ∆t) − c(x, t) ∆t
We can use these approximations to define the numerical approximation to c(x, t). Call this approximation γ, where γm,n represents the approxima- tion evaluated at the point (m∆x, n∆t),
(∂xxc)(m∆x, t) ≈
γm+1,n − 2 γm,n + γm− 1 ,n (∆x)^2
(∂tc) (x, t) ≈
γm,n+1 − γm,n ∆t The mesh approximation of the PDE is thus,
γm,n+1 − γm,n =
D∆t (∆x)^2
{γm+1,n − 2 γm,n + γm− 1 ,n}
3.3 Wave solutions of the mesh equation
Previously, we assumed the exact solution of the diffusion equation took the form of a wave solution. To solve the mesh equation, we assume a solution of the same form, taking m∆x for x and n∆t for t,
γm,n = eσn∆t+ıkm∆x
Evaluating the quantities in the mesh equation,
γm,n+1 − γm,n =
eσ(n+1)∆t^ − eσn∆t
eıkm∆x
eσ∆t^ − 1
γm,n D∆t (∆x)^2
{γm+1,n − 2 γm,n + γm− 1 ,n} =
D∆t (∆x)^2
eık∆x^ − 2 + e−ık∆x
γm,n
The exponentials in parentheses can be simplified using trig identities,
D∆t (∆x)^2
{γm+1,n − 2 γm,n + γm− 1 ,n} = −
4 D∆t (∆x)^2
sin^2
k∆x 2
Comparing these two expressions, and cancelling the factor of γm,n,
eσ∆t^ − 1 =
− 4 D∆t (∆x)^2
sin^2
k∆x 2
This gives the discreet dispersion relation implicitly relating σ to k in the case of the approximate solution.
3.4 Limiting case of the discrete dispersion relation
As a measure of the accuracy of this approximate solution, we expect that this discrete dispersion relation should approximate the exact dispersion relation in the case that the mesh spacing approaches 0. In particular, consider the formal limit as ∆t, ∆x → 0, holding k fixed. Expand the exponential and the sine function using Taylor series,
eσ∆t^ = 1 + σ∆t + O
(∆t)^2
⇒ eσ∆t^ − 1 = σ∆t + O
(∆t)^2
sin
k∆x 2
k∆x 2
+ O
(∆x)^3
⇒ sin^2
k∆x 2
k^2 (∆x)^2 4
+ O
(∆x)^4
4 D∆t (∆x)^2
sin^2
k∆x 2
= −Dk^2 ∆t + O
∆t(∆x)^2
Matching these expressions, the limit of the discrete dispersion relation be- comes,
σ = −Dk^2 + O(∆t) + O
(∆x)^2
→ −Dk^2
Hence the discrete dispersion relation approaches the exact dispersion rela- tion as the mesh size is made very small, as we expect.
Matching real and imaginary components gives two equations,
eS^ cos W = μ eS^ sin W = 0
Since eS^6 = 0, it must be that,
sin W = 0 ⇔ W = nπ
With this value of W , assuming μ < 0, the top equation then becomes,
cos W = (−1)n ⇒ eS^ cos W = eS^ (−1)N^ = −|μ|
Taking N = 1, eS^ = |μ| determines the value of S, which in turn determines σ∆t,
σ∆t = ln |μ| + ıπ
Using the expressions σ∆t = S + ıπ and k = (^) ∆πx , the solution is thus
γm,n = en(σ∆t)+ık(m∆x)^ → enS^ (−1)m+n
If S > 0, then any temporal or spatial mesh oscillations will grow exponen- tially, which is undesirable.
3.5.4 14 < (^) (∆D∆x)t 2 < 12 ⇔ − 1 < μ < 0
In this case, |μ| < 1 ⇒ S < 0. The negative S damps spatial and tempo- ral oscillations, leading to a stable solution. However, due to the complex component of σ∆t, the damping is oscillatory, so the convergence is not as desirable as in the case of (^) (∆D∆x)t 2 < 14.
3.5.5 Summary
In summary, consider the plot of (^) (∆D∆x)t 2 on a number line. Between 0 and 1 4 , there is monotonic decay, which is desirable.^ Between^
1 4 and^
1 2 , there is oscillatory decay, which is acceptable. Greater than 12 , there is oscillatory growth, which is no good.
