Information Technology Reference
In-Depth Information
C
3
above), we can compute the ratio E=x
2
for the four difference cases. The ratio
becomes 0.0633159, 0.0647839, 0.0655636, and 0.0659656. Since these numbers
are almost equal, we conclude that E=x
2
is constant, as it should be. This provides
solid evidence that the program works correctly when the boundary conditions are
u
D
0 and there is no source term. Finding a solution when the boundary conditions
and the source term are non-trivial can be done as described next.
The Method of Manufactured Solutions
We can, in fact, always find an exact solution to a PDE problem. The procedure,
often called the method of manufactured solutions, goes as follows. Pick any func-
tion
v
.x; t / as the exact solution. Insert
v
in the PDE; it does not fit, but adding a
source term will make it fit. As an example, take
v
.x; t /
D
t
2
C
4x
3
:
Inserting this function in the PDE
@t
D
@
2
u
@
u
C
f.x;t/;
@x
2
shows that
v
.x; t / is a solution if we choose
f.x;t/
D
2t
24x:
Then we set I.x/
D
v
.x; 0/, D
0
.t /
D
v
.0; t /,andD
1
.t /
D
v
.1; t /.Totesta
program, we set I , D
0
,andD
1
as just listed and compute E for a sequence of x
values. An approximately constant ratio E=x
2
(for the different E and x values)
provides evidence that the program works.
A nice manufactured solution to use in the debugging phase is a linear function
of space and time. Such a function will very often be (almost) exactly reproduced
by numerical schemes for PDEs, regardless of the values of the grid resolution.
That is, with such a solution we can control that E is sufficiently close to zero (to
machine precision) also on very coarse grids. In the present problem we can choose
v
.x; t /
D
t=2
C
3x as a manufactured solution. The corresponding source term is
f.x;t/
D
1=2. Writing out the error shows that it is about 10
6
,evenoncoarse
grids.
Diffusion of a Jump
Our next problem concerns a jump in the initial condition, like we get when bringing
two metal pieces together, as depicted in Fig.
7.4
. However, for simplicity we apply
Dirichlet conditions (instead of the more appropriate Neumann conditions) at the
boundaries. The problem to be solved reads
