Math578 - Alexiades
1D Diffusion - Comparison with exact solution
Consider a (1-dimensional) diffusion process, with (constant)
diffusivity D > 0, for the concentration u(x,t) in the
semi-infinite interval 0 ≤ x < ∞, during some time 0 ≤ t ≤ tend.
Initialy, the concentration is uniform: uinit(x)=0,
and the endpoint at x=0 is raised and held at concentration u0(t)=1.
This problem admits the explicit similarity solution:
u(x,t) = 1 − ERF( x / (2 SQRT(D t)) ) ,
with ERF(z) the error function.
To test your diffusion code against this exact solution on an interval [0,b],
need to impose the exact solution at x=b:
U(M+1) = 1.0−ERF( b / (2*SQRT(D*t)) ) == ERFC( b / (2*SQRT(D*t)) ) .
So here both BCs are Dirichlet type, and the one at x=b is time-dependent...
think when time needs to be updated...
Take b=1.0, D=0.1, factor=0.96, and compare the numerical and exact solutions
at tend=1, on a mesh with Δx=1/16 (MM=16) and a finer mesh with Δx=1/64 (MM=64).
The comparison should be done in a subprogram COMPARISON, with
minimal changes to the code (so that it can be turned off easily, by commenting out,
on problems without exact solution)
It should print out the numerical solution U(i) at tend, the exact solution u(xi, tend),
and the error (their difference), i.e. print:
xi Ui uEXACT error
Also compute the worst absolute error ( maximum |error| over all i ).
You can plot both U and u_exact on the same plot in gnuplot:
gnuplot>: plot "out" u 1:2 lt 3 lw 2 , "out" u 1:3 lt 1 lw 2
(u 1:2 means using columns 1 and 2 as x and y respectively).
Submit via email as "lab3" :
max errors for MM=16 and MM=64
and the (annotated) plots on paper (or as PDF attachment "lab3_plots.pdf").
CAUTION for fortran programmers: The double precision error function in most fortran
libraries is "DERF", and ERFC is DERFC.
It should be coded as
DERFC( 1.d0 / (2*SQRT(D*t)) )