# Error propagation with NDSolve and periodic boundary conditions

I'm trying to solve the problem of heat conduction through soil with periodic temperature variations on the surface. For now I'm working with a 1D model but I see that the solution doesn't work, it seems that there is some kind of error due to the numerical algorithms used by Mathematica, do you know how to fix that? My code is like that:

a = 20 (*a= years *)
h = 20(*  h= bottom of the simulation*)
s = 365*24*60*60(*s=time in seconds of one year *)
c = 37(* c=cave depth in meters *)
k = 0.756*10^(-6)(*  k=thermal diffusivity*)
d = 20(* d= difference outside and inside temperatures*)
p = (s)(* p= period in seconds of the oscillation*)
delta = 60*60*24*1(*delta= interval for time evaluation*)

sol2 =
NDSolve[
{D[u[t, x], t] == k*D[u[t, x], x, x],
u[0, x] == 0, u[t, 0] == d*Sin[ 2 Pi/p t], D[u[t, h], x] == -u[t, h]},
u, {t, 0, a*s}, {x, 0, h},
Method ->
{"MethodOfLines",
"SpatialDiscretization" ->
{"TensorProductGrid", "DifferenceOrder" -> Pseudospectral"}}];

Print[
"Period of the output temperature oscillation in days:    ",
FunctionPeriod[Sin[2 Pi/p t], t]/60/60/24 // N]

Manipulate[
Plot[u[t, x] /. sol2, {x, 0, 20},
PlotRange -> {-d, d},
ColorFunctionScaling -> False,
ColorFunction -> "TemperatureMap"],
{t, 0, a*s, delta}]

Manipulate[
Plot[u[t, x] /. sol2, {t, 0, a*s},
PlotRange -> Automatic,
ColorFunctionScaling -> False,
ColorFunction -> "TemperatureMap"],
{x, 0, h}]

Manipulate[
Plot[u[t, x] /. sol2, {t, 0, a*s},
PlotRange -> Automatic,
ColorFunctionScaling -> False,
ColorFunction -> "TemperatureMap"],
{{x, 10.65}, 0, 20}]


How you can see from the graph if you move along the x axes the values get wronger and wronger, since I would expect a mean value around 0, instead of this strange behaviour with increased temperature at the beginning. I suspect that the problem is linked with the absolute value of the k constant, in fact setting k=1 we don't have the initial anomaly.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Michael E2 Aug 11 '16 at 13:43
• Thank you for your advice, I took the tour and i found it very useful. – Django Aug 11 '16 at 14:14
• What about these results is wrong? Your graphs show that in the long-term the average temperature at each location seems to be zero, which makes sense given your forcing. The initial dynamics will depend on your initial temperature distribution, which you set constant but that doesn't match the long-term distribution at the beginning of a period. – Chris K Aug 11 '16 at 14:23
• BTW, I wouldn't call these "periodic boundary conditions", which usually implies periodic w/ space so that u[t,0]=u[t,h] and u'[t,0]=u'[t,h]. – Chris K Aug 11 '16 at 14:24
• My intuition is that this solution is in fact correct. Try plotting the spatial distribution of temperature over the first couple cycles with Do[Print[Plot[u[f s, x] /. sol2, {x, 0, h}, PlotRange -> {-20, 20}, PlotLabel -> f]], {f, 0, 2, 0.1}]; In the first cycle, the surface initially warms and this temperature diffuses down, so the deeper temperature has no choice but to increase. After a few cycles, you see that the temperature distribution at the beginning of a period is not constant as in your initial conditions. – Chris K Aug 11 '16 at 14:57