The 1D Diffusion Equation
The 1D Diffusion Equation
The 1D Diffusion Equation
=
t
where u(x, t) is the unknown function to be solved for, x is a coordinate in space, and t is time. The
coefficient isthediffusioncoefficientanddetermineshowfastuchangesintime.Aquickshortformforthe
diffusionequationisut = uxx .
Comparedtothewaveequation, utt = c2 uxx ,whichlooksverysimilar,butthediffusionequationfeatures
solutions that are very different from those of the wave equation. Also, the diffusion equation makes quite
differentdemandstothenumericalmethods.
Typical diffusion problems may experience rapid change in the very beginning, but then the evolution of u
becomesslowerandslower.Thesolutionisusuallyverysmooth,andaftersometime,onecannotrecognize
the initial shape of u. This is in sharp contrast to solutions of the wave equation where the initial shape is
2
u
(x) = 0 .ThisstationarylimitofthediffusionequationiscalledtheLaplaceequationandarisesinavery
widerangeofapplicationsthroughoutthesciences.
Itispossibletosolvefor u(x, t) usingaexplicitscheme,butthetimesteprestrictionssoonbecomemuch
lessfavorablethanforanexplicitschemeforthewaveequation.Andofmoreimportance,sincethesolutionu
of the diffusion equation is very smooth and changes slowly, small time steps are not convenient and not
requiredbyaccuracyasthediffusionprocessconvergestoastationarystate.
(1)
=
t
x [0, L]
(2)
u(0, t) = 0,
t > 0,
(3)
u(L, t) = 0,
t > 0.
(4)
u(x, 0) = I(x),
Equation (28) is known as a onedimensional diffusion equation, also often referred to as a heat equation.
Withonlyafirstorderderivativeintime,onlyoneinitialconditionisneeded,whilethesecondorderderivative
inspaceleadstoademandfortwoboundaryconditions.Theparameter mustbegivenandisreferredtoas
thediffusioncoefficient.
Diffusionequationslike(28) have a wide range of applications throughout physical, biological, and financial
sciences. One of the most common applications is propagation of heat, where u(x, t) represents the
temperatureofsomesubstanceatpointxandtimet .Thesectiondiffu:appgoesintoseveralwidelyoccurring
applications.
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
1/23
10/02/2015
[0, T ] byasetofmeshpoints.
i = 0, , N x ,
and
t n = nt,
Moreover,
n
i
n = 0, , N t .
and
n = 0, , N t
u(xi , t n ) =
(5)
u(xi , t n ),
Thenextstepistoreplacethederivativesbyfinitedifferenceapproximations.Thecomputationallysimplest
methodarisesfromusingaforwarddifferenceintimeandacentraldifferenceinspace:
+
[D
u = Dx Dx u]
n
i
(6)
Writtenout,
u
n+1
i
n
i
u
=
n
i+1
2u
+ u
(7)
n
i1
WehaveturnedthePDEintoalgebraicequations,alsooftencalleddiscreteequations.Thekeypropertyof
theequationsisthattheyarealgebraic,whichmakesthemeasytosolve.Asusual,weanticipatethat un
is
i
alreadycomputedsuchthatun+1
istheonlyunknownin(7).Solvingwithrespecttothisunknowniseasy:
i
u
n+1
i
= u
n
i
+ F (u
n
i+1
2u
n
i
+ u
n
i1
) .
(8)
F isthekeyparameterinthediscretediffusionequation:
NotethatF isadimensionlessnumber
thatlumpsthekeyphysicalparameterintheproblem, ,andthediscretizationparametersx andt
intoasingleparameter.Allthepropertiesofthenumericalmethodarecriticallydependentuponthevalueof
F (seethesectiondiffu:pde1:analysisfordetails).
Thecomputationalalgorithmthenbecomes
1. compute$u^0_i=I(x_i)$fori
2. forn = 0, 1, , Nt :
= 0, , N x
1. apply(8)foralltheinternalspatialpointsi
n+1
2. settheboundaryvaluesui
= 1, , N x 1
= 0 fori = 0 andi = N x
ThealgorithmiscompactlyfullyspecifiedinPython:
x = linspace(0, L, Nx+1)
dx = x[1] - x[0]
t = linspace(0, T, Nt+1)
dt = t[1] - t[0]
F = a*dt/dx**2
u = zeros(Nx+1)
u_1 = zeros(Nx+1)
2/23
10/02/2015
Theprogramdiffu1D_u0.pycontainsafunction solver_FEforsolvingthe1Ddiffusionequationwithu = 0 on
theboundary.Thefunctions plugand gaussianrunsthecasewith I(x) asadiscontinuousplugorasmooth
Gaussianfunction,respectively.Experimentswiththesetwofunctionsrevealsomeimportantobservations:
TheForwardEulerschemeleadstogrowingsolutionsifF
>
I(x) asadiscontinuousplugleadstoasawtoothlikenoisefor F =
isabsentforF
1
4
1
2
,seemovie,which
,seemovie.
ThesmoothGaussianinitialfunctionleadstoasmoothsolution,seemovieforF
1
2
[D u = Dx Dx u] ,
t
whichwrittenoutreads
u
n
i
n1
i
u
=
n
i+1
2u
n
i
+ u
(10)
n
i1
Nowweassume un1
iscomputed,butallquantitiesatthenewtimelevel n areunknown.Thistimeitis
i
notpossibletosolvewithrespecttoun
becausethisvaluecouplestoitsneighborsinspace,un
andun
,
i
i1
i+1
which are also unknown. Let us examine this fact for the case when Nx = 3 . Equation (10) written for
i = 1, , N x 1 = 1, 2 becomes
u
n
1
n1
1
u
=
n
2
2u
t
u
n
2
x
n1
2
u
=
n
3
2u
+ u
+ u
n
1
Theboundaryvalues un
and un
areknownaszero.Collectingtheunknownnewvalues un
and un
onthe
0
3
1
2
lefthandsidegives
(1 + 2F ) u
F u
Thisisacoupled 2
formis
n
1
n
1
Fu
n
2
+ (1 + 2F ) u
n
2
= u
n1
1
= u
n1
2
2systemofalgebraicequationsfortheunknowns u
1 + 2F
)(
F
1 + 2F
n
1
n
2
u
) = (
u
n
1
and un
.Theequivalentmatrix
2
n1
1
n1
Implicitvs.explicitmethods: Discretizationmethodsthatleadtoacoupledsystemofequationsforthe
unknownfunctionatanewtimelevelaresaidtobeimplicitmethods.Thecounterpart,explicitmethods,
referstodiscretizationmethodswherethereisasimpleexplicitformulaforthevaluesoftheunknown
functionateachofthespatialmeshpointsatthenewtimelevel.Fromanimplementationalpointofview,
implicitmethodsaremorecomprehensivetocodesincetheyrequirethesolutionofcoupledequations,i.e.,
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
3/23
10/02/2015
amatrixsystem,ateachtimelevel.
Inthegeneralcase,(10)givesrisetoacoupled(N x 1) (N x 1)systemofalgebraicequationsfor
alltheunknown un
attheinteriorspatialpoints i = 1, , N x 1 .Collectingtheunknownsontheleft
i
handside,(10)canbewritten
F u
fori
n
i1
+ (1 + 2F ) u
n
i
Fu
n
i+1
= u
n1
i1
(11)
= 1, , N x 1 .Here,wehaveintroducedthemeshFouriernumber
t
F =
x
One can either view these equations as a system for where the un
values at the internal mesh points,
i
n
n
i = 1, , N x 1 ,areunknown,orwemayappendtheboundaryvalues u and u
tothesystem.In
0
N
the latter case, all un
for i
i
N x 1 equationsin(11):
u
u
n
0
n
Nx
(12)
= 0,
(13)
= 0.
Acoupledsystemofalgebraicequationscanbewrittenonmatrixform,andthisisimportantifwewanttocall
up readymade software for solving the system. Theequations(11)and(12)(13) correspond to the matrix
equation
AU = b
whereU
Nx
= (u , , u
A 0,0
A 1,0
A =
) ,andthematrixA hasthefollowingstructure:
A 0,1
A 1,1
A 2,1
A 2,2
A 2,3
A i,i1
A i,i
A i,i+1
A Nx 1,Nx
A Nx ,Nx 1
A N x ,N x
(14)
Thenonzeroelementsaregivenby
A i,i1 = F
A i,i = 1 + 2F
A i,i+1 = F
fortheequationsforinternalpoints,i
to
= 1, , N x 1 .Theequationsfortheboundarypointscorrespond
A 0,0 = 1,
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
4/23
10/02/2015
Therighthandsidebiswrittenas
b0
b1
b =
bi
bN x
with
b0 = 0,
bi = u
n1
i
i = 1, , N x 1,
bN x = 0 .
We observe that the matrix A contains quantities that do not change in time. Therefore, A can be formed
onceandforallbeforeweentertherecursiveformulasforthetimeevolution.Therighthandsideb,however,
mustbeupdatedateachtimestep.Thisleadstothefollowingcomputationalalgorithm,heresketchedwith
Pythoncode:
x = linspace(0, L, Nx+1)
dx = x[1] - x[0]
t = linspace(0, T, N+1)
u = zeros(Nx+1)
u_1 = zeros(Nx+1)
5/23
10/02/2015
Wehaveseenfrom(14)thatthematrix A istridiagonal.Thecodesegmentaboveusedafull,densematrix
representation of A , which stores a lot of values we know are zero beforehand, and worse, the solution
algorithm computes with all these zeros. With Nx + 1 unknowns, the work by the solution algorithm is
1
and the storage requirements (Nx + 1) . By utilizing the fact that A is tridiagonal and
employingcorrespondingsoftwaretools,theworkandstoragedemandscanbeproportionaltoNx only.
3
(N x + 1)
The key idea is to apply a data structure for a tridiagonal or sparse matrix. The scipy.sparse package has
relevantutilities.Forexample,wecanstorethenonzerodiagonalsofamatrix.Thepackagealsohaslinear
system solvers that operate on sparse matrix data structures. The code below illustrates how we can store
onlythemaindiagonalandtheupperandlowerdiagonals.
# Representation of sparse matrix and right-hand side
main = zeros(Nx+1)
lower = zeros(Nx-1)
upper = zeros(Nx-1)
b
= zeros(Nx+1)
# Precompute sparse matrix
main[:] = 1 + 2*F
lower[:] = -F #1
upper[:] = -F #1
# Insert boundary conditions
main[0] = 1
main[Nx] = 1
A = scipy.sparse.diags(
diagonals=[main, lower, upper],
offsets=[0, -1, 1], shape=(Nx+1, Nx+1),
format='csr')
print A.todense() # Check that A is correct
# Set initial condition
for i in range(0,Nx+1):
u_1[i] = I(x[i])
for n in range(0, Nt):
b = u_1
b[0] = b[-1] = 0.0 # boundary conditions
u[:] = scipy.sparse.linalg.spsolve(A, b)
u_1[:] = u
Crank-Nicolson scheme
TheideaintheCrankNicolsonschemeistoapplycentereddifferencesinspaceandtime,combinedwithan
averageintime.WedemandthePDEtobefulfilledatthespatialmeshpoints,butinbetweenthepointsinthe
timemesh:
fori
u(xi , t
n+
) =
u(xi , t
n+
).
= 1, , N x 1 andn = 0, , N t 1 .
Withcentereddifferencesinspaceandtime,weget
n+
[Dt u = Dx Dx u]
1
2
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
6/23
10/02/2015
Ontherighthandsidewegetanexpression
1
n+
(u
n+
Thisexpressionisproblematicsinceui
n+
n+
2u
i1
n+
+ u
1
2
i+1
) .
1
2
isnotoneoftheunknownwecompute.Apossibilityistoreplace
1
2
byanarithmeticaverage:
n+
(u
n
i
+ u
n+1
i
) .
Inthecompactnotation,wecanusethearithmeticaveragenotation
u
:
[Dt u = Dx Dx
u
]
n+
1
2
Afterwritingoutthedifferencesandaverage,multiplyingby t ,andcollectingallunknowntermsontheleft
handside,weget
u
n+1
i
F (u
2
n+1
i1
2u
n+1
i
+ u
n+1
i+1
) = u
n
i
1
+
F (u
2
n
i1
2u
n
i
+ u
n
i+1
).
Alsohere,asintheBackwardEulerscheme,thenewunknowns un+1
, un+1
,and un+1
arecoupledina
i1
i
i+1
linearsystemAU
= b,whereA hasthesamestructureasin(14),butwithslightlydifferententries:
1
A i,i1 =
F
2
1
A i,i =
+ F
2
1
A i,i+1 =
fortheequationsforinternalpoints,i
to
F
2
= 1, , N x 1 .Theequationsfortheboundarypointscorrespond
A 0,0 = 1,
A 0,1 = 0,
A Nx ,Nx 1 = 0,
A N x ,N x = 1 .
Therighthandsidebhasentries
b0 = 0,
bi = u
n1
i
i = 1, , N x 1,
bN x = 0 .
The
rule
Theruleprovidesafamilyoffinitedifferenceapproximationsintime:
= 0 givestheForwardEulerschemeintime
= 1 givestheBackwardEulerschemeintime
1
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
=
7/23
10/02/2015
1
2
givestheCrankNicolsonschemeintime
Appliedtothe1Ddiffusionproblemwehave
u
n+1
i
n
i
u
= (
n+1
2u
i+1
n+1
i
+ u
n+1
i1
u
+ (1 )
n
i+1
2u
n
i
+ u
n
i1
) .
Thisschemealsoleadstoamatrixsystemwithentries
A i,i1 = F ,
A i,i = 1 + 2F
, A i,i+1 = F ,
whilerighthandsideentrybi is
bi = u
n
i
u
+ F (1 )
n
i+1
2u
n
i
+ u
n
i1
ThecorrespondingentriesfortheboundarypointsareasintheBackwardEulerandCrankNicolsonschemes
listedearlier.
throughout science and engineering. In 1D these equations read u (x) = 0 and u (x) = f(x) ,
respectively. We can solve 1D variants of the Laplace equations with the listed software, because we can
interpretuxx = 0 asthelimitingsolutionof ut = uxx when ureachasteadystatelimitwhere ut 0 .
Similarly,Poissonsequationuxx = f arisesfromsolvingut = uxx + f andlettingt sout 0 .
Technicallyinaprogram,wecansimulate t byjusttakingonelargetimestep,orequivalently,set
toalargevalue.AllweneedistohaveF large.AsF ,wecanfromtheschemesseethatthelimiting
discreteequationbecomes
u
n+1
i+1
2u
n+1
i
whichisnothingbutthediscretization[Dx Dx u]n+1
i
+ u
n+1
i1
= 0,
= 0 ofuxx = 0 .
The Backward Euler scheme can solve the limit equation directly and hence produce a solution of the 1D
Laplace equation. With the Forward Euler scheme we must do the time stepping since F > 1/2 is illegal
and leads to instability. We may interpret this time stepping as solving the equation system from uxx by
iteratingonatimepseudotimevariable.
Extensions
These extensions are performed exactly as for a wave equation as they only affect the spatial derivatives
(whicharethesameasinthewaveequation).
Variablecoefficients
NeumannandRobinconditions
2Dand3D
Future versions of this document will for completeness and independence of the wave equation document
featureinfoonthethreepoints.TheRobinconditionisnew,butstraightforwardtohandle:
u
= hT (u Us ),
[Dx u = hT (u Us )]
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
n
i
8/23
10/02/2015
ut = uxx ,
is that the initial shape u(x, 0) = I(x) spreads out in space with time, along with a decaying amplitude.
Threedifferentexampleswillillustratethespreadingofuinspaceandthedecayintime.
Similarity solution
The diffusion equation (15) admits solutions that depend on
Oneparticularsolutionis
(16)
u(x, t) = a erf() + b,
where
erf() =
(17)
d,
erf() = 1,
lim erf() = 1,
erf() = erf(),
erf(0) = 0,
erf(2) = 0.99532227,
erf(3) = 0.99997791 .
As t 0 ,theerrorfunctionapproachesastepfunctioncenteredat x = c .Foradiffusionproblemposed
on the unit interval [0, 1] , we may choose the step at x = 1/2 (meaning c = 1/2 ), a = 1/2,
b = 1/2 .Then
x
1
(1 erf (
u(x, t) =
2
(18)
1
2
erfc (
)) =
),
2
4t
4t
= 1 erf() .
1/2
1
(1 erf (
u(0, t) =
1/2
1
(1 erf (
u(1, t) =
2
(19)
)) ,
4t
(20)
))
4t
Forsmallenought ,u(0, t)
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
9/23
10/02/2015
= uxx admitsaGaussianfunctionassolution:
u(x, t) =
(21)
(x c)
exp (
4t
)
4t
At t = 0 thisisaDiracdeltafunction,soforcomputationalpurposesonemuststarttoviewthesolutionat
sometime t = t > 0 .Replacing t by t + t in(21)makesiteasytooperatewitha(new) t thatstartsat
t = 0 withaninitialconditionwithafinitewidth.Theimportantfeatureof(21)isthatthestandarddeviation
of a sharp initial Gaussian pulse increases in time according to
flattenout.
at
(22)
sin (kx)
TheparametersQandk canbefreelychosen,whileinserting(22)in(15)givestheconstraint
a = k
A very important feature is that the initial shape I(x) = Q sin kx undergoes a damping exp (k2 t) ,
meaningthatrapidoscillationsinspace,correspondingtolarge k ,areverymuchfasterdampenedthanslow
oscillationsinspace,correspondingtosmall k .Thisfeatureleadstoasmoothingoftheinitialconditionwith
time.
Thefollowingexamplesillustratesthedampingpropertiesof(22).Weconsiderthespecificproblem
ut = uxx ,
u(0, t) = u(1, t) = 0,
t (0, T ],
Theinitialconditionhasbeenchosensuchthataddingtwosolutionslike(22)constructsananalyticalsolution
totheproblem:
u(x, t) = e
sin(x) + 0.1e
10 t
sin(100x)
(23)
FigureEvolutionofthesolutionofadiffusionproblem:initialcondition(upperleft),1/100reductionofthesmall
waves(upperright),1/10reductionofthelongwave(lowerleft),and1/100reductionofthelongwave(lower
right)illustratestherapiddampingofrapidoscillationssin(100x)andtheverymuchslowerdampingofthe
4
slowly varying sin(x) term. After about t = 0.5 10 the rapid oscillations do not have a visible
amplitude,whilewehavetowaituntilt 0.5 beforetheamplitudeofthelongwavesin(x) becomesvery
small.
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
10/23
10/02/2015
Evolutionofthesolutionofadiffusionproblem:initialcondition(upperleft),1/100reductionofthesmallwaves
(upperright),1/10reductionofthelongwave(lowerleft),and1/100reductionofthelongwave(lowerright)
UL ,
x < L/2
UR ,
x L/2
Suchadiscontinuousinitialconditionmayarisewhentwoinsulatedblocksofmetalsatdifferenttemperature
are brought in contact at t = 0 . Alternatively, signaling in the brain is based on release of a huge ion
concentration on one side of a synapse, which implies diffusive transport of a discontinuous concentration
function.
Moretobewritten...
at
ikx
wherei = 1 istheimaginaryunit.Wecanaddsuchfunctions,oftenreferredtoaswavecomponents,to
makeaFourierrepresentationofageneralsolutionofthediffusionequation:
t ikx
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
u(x, t)
,
k
11/23
10/02/2015
u(x, t) bk e
k t
ikx
(24)
kK
ikx
= 0 we
kK
(The relevant formulas for bk come from Fourier analysis, or equivalently, a leastsquares method for
approximatingI(x) inafunctionspacewithbasisexp (ikx) .)
Much insight about the behavior of numerical methods can be obtained by investigating how a wave
component
components are also solutions of the schemes, but the damping factor exp (k t) varies among the
schemes.Toeasetheforthcomingalgebra,wewritethedampingfactorasAn .Theexactamplificationfactor
correspondingtoA isAe
= exp (k t) .
k t
ikx
Afundamentalquestioniswhethersuchcomponentsarealsosolutionsofthefinitedifferenceschemes.This
2
isindeedthecase,buttheamplitudeexp (k t) mightbemodified(whichalsohappenswhensolvingthe
ODEcounterpartu = u ).Wethereforelookfornumericalsolutionsoftheform
n
uq = A
ikqx
= A
ikx
(25)
wheretheamplificationfactorA mustbedeterminedbyinsertingthecomponentintoanactualscheme.
Stability (1)
2
Accuracy (1)
To determine how accurately a finite difference scheme treats one wave component (25), we see that the
n
basic deviation from the exact solution is reflected in how well An approximates Ae , or how well A
approximatesAe .
[D
= uxx canbewrittenas
n
u = Dx Dx u]q .
Insertingawavecomponent(25)intheschemedemandscalculatingtheterms
A 1
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
ikqx
n
ikqx
n
[
12/23
10/02/2015
ikqx
[D
A]
= e
ikqx
A 1
,
t
and
A
Dx Dx [e
ikx
]q = A
(e
ikqx
sin
kx
(
)) .
2
4
=
sin
kx
(
),
2
andconsequently
A = 1 4F sin
kx
(
),
2
where
t
F =
x
isthenumericalFouriernumber.Thecompletenumericalsolutionisthen
n
uq
= (1 4F sin
kx
(
))
ikqx
Stability (2)
Weeasilyseethat A 1 .However,the A canbelessthan 1,whichwillleadtogrowthofanumerical
wavecomponent.ThecriterionA 1 implies
2
4F sin (p/2) 2 .
Theworstcaseiswhensin2 (p/2)
= 1 ,soasufficientcriterionforstabilityis
1
F
,
2
orexpressedasaconditionont :
x
.
2
Accuracy (2)
SinceA isexpressedintermsofF andtheparameterwenowcallp
andp :
= kx/2 ,wealsoexpressA e by F
13/23
10/02/2015
Theresultis
A
= 1 4F sin
p + 2F p
16F
sin
p + 8F
Ae
Recalling that F =
termsareatmost
t
1 4
x
+ t 4 t
+ t x
+ .
= uxx byaBackwardEulerscheme,
n
[D u = Dx Dx u]q ,
t
and inserting a wave component (25), leads to calculations similar to those arising from the Forward Euler
scheme,butsince
ikqx
[D A]
t
= A
ikqx
1 A
,
t
weget
1 A
4
=
sin
kx
),
2
andthen
A = (1 + 4F sin
(26)
p)
Thecompletenumericalsolutioncanbewritten
n
uq = (1 + 4F sin
p)
ikqx
Stability (3)
Weseefrom(26)that 0
oscillatoryforanyt >
n+
[Dt u = Dx Dx
u
]q
1
2
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
14/23
10/02/2015
or
n+
[Dt u]q
=
2
n+1
) .
Inserting(25)inthetimederivativeapproximationleadsto
1
[Dt A
ikqx n+
1
2
= A
n+
1
2
ikqx
1
2
= A
ikqx
A 1
.
t
1
=
sin
kx
(
) (1 + A),
2
andaftersomemorealgebra,
1 2F sin
A =
1 + 2F sin
p
.
p
Theexactnumericalsolutionishence
n
uq
1 2F sin
= (
1 + 2F sin
p
)
ikpx
Stability (4)
ThecriteriaA
Amplificationfactorsforlargetimesteps
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
15/23
10/02/2015
AmplificationfactorsfortimestepsaroundtheForwardEulerstabilitylimit
Amplificationfactorsforsmalltimesteps
Theeffectofnegativeamplificationfactorsisthat An changessignfromonetimeleveltothenext,thereby
givingrisetooscillationsintimeinananimationofthesolution.WeseefromFigureAmplificationfactorsfor
largetimestepsthatfor F = 20 ,waveswith p /2 undergoadampingcloseto 1,whichmeansthat
theamplitudedoesnotdecayandthatthewavecomponentjumpsupanddownintime.ForF = 2 wehave
a damping of a factor of 0.5 from one time level to the next, which is very much smaller than the exact
damping.Shortwaveswillthereforefailtobeeffectivelydampened.Thesewaveswillmanifestthemselvesas
highfrequencyoscillatorynoiseinthesolution.
ikx
=. . . ...
16/23
10/02/2015
L
2
e dx
0
E(t) =
.
L
udx
0
Isthisasuitableerrormeasureforthepresentproblem?
d)Insteadofusing u(L, t) = 0 asapproximateboundaryconditionforlettingthediffusedGaussianpulse
outofourfinitedomain,onemaytry ux (L, t) = 0 sincethesolutionforlarge t isquiteflat.Arguethatthis
conditiongivesacompletelywrongasymptoticsolutionast 0 .Todothis,integratethediffusionequation
from 0 to L , integrate uxx by parts (or use Gauss divergence theorem in 1D) to arrive at the important
property
L
dt
L
implyingthat0
u(x, t)dx = 0,
0
udx mustbeconstantintime,andtherefore
L
u(x, t)dx =
0
I(x)dx .
0
Theintegraloftheinitialpulseis1.
e)Anotherideaforanartificialboundaryconditionatx
= L istouseacoolinglaw
ux = q(u uS ),
(27)
17/23
10/02/2015
diffu_symmetric_gaussian_2D.pdf.
a)DeriveindetailaForwardEulerscheme,aBackwardEulerscheme,andaCrankNicolsonforthistypeof
diffusionmodel.Thereafter,formulatearuletosummarizethethreeschemes.
b)Assumeasolutionlike(22)andfindtherelationbetweena ,k , ,and .
c)CalculatethestabilityoftheForwardEulerscheme.Designnumericalexperimentstoconfirmtheresults.
d)Repeatc)fortheBackwardEulerscheme.
e)Repeatc)fortheCrankNicolsonscheme.
f)Howdoestheextratermbu impacttheaccuracyofthethreeschemes?
Hint.Comparethenumericalandexactamplificationfactors,eitheringraphsorbyTaylorseriesexpansion
(orboth).
Filename: diffu_stab_uterm.pdf.
(28)
u
((x)
= (x). A 1D
),
u(x, 0) = I(x),
(29)
x [0, L]
u(0, t) = U0 ,
u(L, t) = UL ,
t > 0,
(30)
t > 0.
(31)
Ashortformofthediffusionequationwithvariablecoefficientsisut
= (ux )x .
Stationary solution
As t , the solution of the above problem will approach a stationary limit where u/t
governingequationisthen
d
dx
The
(32)
du
(
= 0 .
) = 0,
dx
u(x) = U0 + (UL U0 )
x
0
L
(())
(33)
d
.
(())
18/23
10/02/2015
Consider a medium built of M layers. The boundaries between the layers are denoted by b0 , , bM ,
whereb0 = 0and bM = L .Ifthematerialineachlayerpotentiallydiffersfromtheothers,butisotherwise
constant,wecanexpress asapiecewiseconstantfunctionaccordingto
0 ,
(x) = i ,
b0 x < b1 ,
(34)
bi x < bi+1 ,
0 ,
bM1 x bM .
Theexactsolution(33)incaseofsuchapiecewiseconstant functioniseasytoderive.Assumethatxisin
x
m1
j=0
d wemustintegratethroughthefirst m 1
bm intothem thlayer:
u(x) = U0 + (UL U0 )
M1
j=0
(35)
(bj+1 bj )/(bj )
Remark. It may sound strange to have a discontinuous in a differential equation where one is to
differentiate,butadiscontinuous iscompensatedbyadiscontinuous ux suchthat ux iscontinuesand
thereforecanbedifferentiatedas(ux )x .
Implementation
Programmingwithpiecewisefunctiondefinitionquicklybecomescumbersomeasthemostnaiveapproachis
totestforwhichintervalxlies,andthenstartevaluatingaformulalike(35).InPython,vectorizedexpressions
may help to speed up the computations. The convenience classes PiecewiseConstant and
IntegratedPiecewiseConstantweremadetosimplifyprogrammingwithfunctionslike(34)andexpressionslike
(35). These utilities not only represent piecewise constant functions, but also smoothed versions of them
where the discontinuities can be smoothed out in a controlled fashion. This is advantageous in many
computationalcontexts(althoughseldomforpurefinitedifferencecomputationsofthesolutionu).
The PiecewiseConstant class is created by sending in the domain as a 2tuple or 2list and a data object
describing the boundaries b0 , , bM and the corresponding function values 0 , , M1 . More
precisely, dataisanestedlist,where data[i][0]holds bi and data[i][1]holdsthecorrespondingvalue i ,
fori = 0, , M 1 .Givenbi andi inarrays band a,itiseasytofilloutthenestedlist data.
Inourapplication,wewanttorepresent and 1/ aspiecewiseconstantfunction,inadditiontothe u(x)
function which involves the integrals of 1/ . A class creating the functions we need and a method for
evaluatingu,cantaketheform
class SerialLayers:
"""
b: coordinates of boundaries of layers, b[0] is left boundary
and b[-1] is right boundary of the domain [0,L].
a: values of the functions in each layer (len(a) = len(b)-1).
U_0: u(x) value at left boundary x=0=b[0].
U_L: u(x) value at right boundary x=L=b[0].
"""
def __init__(self, a, b, U_0, U_L, eps=0):
self.a, self.b = np.asarray(a), np.asarray(b)
self.eps = eps # smoothing parameter for smoothed a
self.U_0, self.U_L = U_0, U_L
a_data = [[bi, ai] for bi, ai in zip(self.b, self.a)]
domain = [b[0], b[-1]]
self.a_func = PiecewiseConstant(domain, a_data, eps)
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
19/23
10/02/2015
Avisualizationmethodisalsoconvenienttohave.Belowweplotu(x)alongwith(x) (whichworkswellas
longasmax (x) isofthesamesizeasmax u = max(U0 , UL ) ).
class SerialLayers:
...
def plot(self):
x, y_a = self.a_func.plot()
x = np.asarray(x); y_a = np.asarray(y_a)
y_u = self.u_exact(x)
import matplotlib.pyplot as plt
plt.figure()
plt.plot(x, y_u, 'b')
plt.hold('on') # Matlab style
plt.plot(x, y_a, 'r')
ymin = -0.1
ymax = 1.2*max(y_u.max(), y_a.max())
plt.axis([x[0], x[-1], ymin, ymax])
plt.legend(['solution $u$', 'coefficient $a$'], loc='upper left')
if self.eps > 0:
plt.title('Smoothing eps: %s' % self.eps)
plt.savefig('tmp.pdf')
plt.savefig('tmp.png')
plt.show()
FigureSolutionofthestationarydiffusionequationcorrespondingtoapiecewiseconstantdiffusioncoefficient
showsthecasewhere
b = [0, 0.25, 0.5, 1]
a = [0.2, 0.4, 4]
U_0 = 0.5; U_L = 5
# material boundaries
# material values
# boundary conditions
Solutionofthestationarydiffusionequationcorrespondingtoapiecewiseconstantdiffusioncoefficient
Byaddingthe epsparametertotheconstructorofthe SerialLayersclass,wecanexperimentwithsmoothed
versions of and see the (small) impact on u. Figure Solution of the stationary diffusion equation
correspondingtoa*smoothedpiecewiseconstantdiffusioncoefficient*showstheresult.
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
20/23
10/02/2015
Solutionofthestationarydiffusionequationcorrespondingtoa*smoothedpiecewiseconstantdiffusion
coefficient*
Exercises
Exercise 4: Stabilizing the Crank-Nicolson method by Rannacher time
stepping
It is well known that the CrankNicolson method may give rise to nonphysical oscillations in the solution of
diffusionequationsiftheinitialdataexhibitjumps(seethesectionAnalysis of the CrankNicolson scheme).
Rannacher[Ref1] suggested a stabilizing technique consisting of using the Backward Euler scheme for the
firsttwotimestepswithsteplength 1 t.Onecangeneralizethisideatotaking2m timestepsofsize 1 t
2
2
with the Backward Euler method and then continuing with the CrankNicolson method, which is of second
order in time. The idea is that the high frequencies of the initial solution are quickly damped out, and the
BackwardEulerschemetreatsthesehighfrequenciescorrectly.Thereafter,thehighfrequencycontentofthe
solutionisgoneandtheCrankNicolsonmethodwilldowell.
Test this idea for m = 1, 2, 3 on a diffusion problem with a discontinuous initial condition. Measure the
convergence rate using the solution (18) with the boundary conditions (19)(20) for t values such that the
conditionsareinthevicinityof1.Forexample,t < 5a1.6 102 makesthesolutiondiffusionfromastep
toalmostastraightline.Theprogram diffu_erf_sol.pyshowshowtocomputetheanalyticalsolution.
(36)
(37)
(38)
Theenergyestimateforthisproblemreads
||u||
wherethe||
||
||I||
(39)
normisdefinedby
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
||g|
21/23
10/02/2015
(40)
||g||
g dx .
0
Theestimate(39)saysthatthesizeof$u$neverexceedsthatoftheinitialcondition,ormoreequivalently,
thattheareaundertheucurvedecreaseswithtime.
Toshow(39),multiplythePDEby uandintegratefrom 0 to L .Usethat uut canbeexpressedasthetime
derivative of u2 and that ux xu can integrated by parts to form an integrand u2x . Show that the time
2
derivativeof||u||L mustbelessthanorequaltozero.Integratethisexpressionandderive(39).
2
b)Nowweaddressaslightlydifferentproblem,
ut = ux x + f(x, t), x = (0, L), t (0, T ],
(41)
(42)
u(x, 0) = 0, x [0, L]
(43)
Theassociatedenergyestimateis
||u||
||f||
(44)
(Thisresultismoredifficulttoderive.)
NowconsiderthecompoundproblemwithaninitialconditionI(x) andarighthandsidef(x, t) :
ut = ux x + f(x, t), x = (0, L), t (0, T ],
(45)
(46)
(47)
showthattheenergyestimatefor(45)(47)becomes
||u||
||I||
+ ||f||
(48)
c)Oneapplicationof(48) is to prove uniqueness of the solution. Suppose u1 and u2 both fulfill (45)(47).
Showthat u = u1 u2 thenfulfills(45)(47)with f = 0and I = 0.Use(48)todeducethattheenergy
mustbezeroforalltimesandthereforethatu1 = u2 ,whichprovesthatthesolutionisunique.
d)Generalize(48)toa2D/3Ddiffusionequationut
= (u) forx .
Hint.Useintegrationbypartsinmultidimensions:
u
u (u) dx =
where n
u u dx +
,
n
= n u ,n beingtheoutwardunitnormaltotheboundary ofthedomain .
e)NowwealsoconsiderthemultidimensionalPDEut
q dx =
q n ds
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
22/23
10/02/2015
q n ds
= 0 ,
u dx =
I dx .
f)Establishacodein1D,2D,or3Dthatcansolveadiffusionequationwithasourceterm f ,initialcondition
I ,andzeroDirichletorNeumannconditionsonthewholeboundary.
Wecanuse(48)and(49)asapartialverificationofthecode.Choosesomefunctionsf andI andcheckthat
(48) is obeyed at any time when zero Dirichlet conditions are used. Iterate over the same I functions and
checkthat(49)isfulfilledwhenusingzeroNeumannconditions.
g)Makealistofsomepossiblebugsinthecode,suchasindexingerrorsinarrays,failuretosetthecorrect
boundaryconditions,evaluationofatermatawrongtimelevel,andsimilar.Foreachofthebugs,seeifthe
verificationtestsfromtheprevioussubexercisepassorfail.Thisinvestigationshowshowstrongtheenergy
estimatesandtheestimate(49)areforpointingouterrorsintheimplementation.
Filename: diffu_energy.pdf.
http://hplgit.github.io/num-methods-for-PDEs/doc/pub/diffu/sphinx/._main_diffu001.html
23/23