# Solution to heat equation with Neumann boundary conditions and lateral heat loss

I am conducting an experiment on heat conduction in a copper rod. To do so I am using a thin copper rod with several temperature probes, then I am going to compare the data with a solution to the heat equation accounting for lateral heat loss. The bar is $0.2 \textrm{m}$ long and the thermal diffusivity ($\alpha$) of copper is $1.12*10^{-4} m^2 s^{-1}$ and the coefficient of lateral heat loss ($\kappa$) is $\ 2.61*10^{-3} s^{-1}$. First I calculated the initial temperature distribution of the bar and got:

  data =
{{0, 26.75}, {0.01, 26.5}, {0.02, 28}, {0.03, 28.75}, {0.04, 30}, {0.05,  32},
{0.06, 34.25}, {0.07, 37.5}, {0.08, 43.25}, {0.09, 59.5}, {0.1, 102},
{0.11, 67.75}, {0.12, 40.5}, {0.13, 35.75}, {0.14, 34}, {0.15, 31.75},
{0.16, 30}, {0.17, 27.5}, {0.18, 28}, {0.19, 26.5}, {0.2, 28}}


The graph of this data:

From here, I tried to use Mathematica to find a Fourier series with help from this. And I also calculated it by hand which is:

0.7607201329841713/2 +
Sum[(1/10 ((0.912752572216549 Sin[3.141592653589793 n])/
n + (-0.002189478133985855 +
5.795634757369874 Cos[3.141592653589793 n] +
2.310248484743582 n Sin[3.141592653589793 n])/(
6.293389749792516 +
1. n^2) + (-0.888276236822269 Sin[3.141592653589793 n] +
0.888276236822269 Sin[6.283185307179586 n])/n +
1/(5.3549458751438115 +
1. n^2) (5.51162427937814 Cos[3.141592653589793 n] -
0.003837157071497423 Cos[6.283185307179586 n] -
2.3817822905923016 n Sin[3.141592653589793 n] +
0.0016581813809965903 n Sin[6.283185307179586 n])))*
Cos[10 (n) π (x)] + (1/
10 ((5.74431 +
0.913625 n^2 + (-5.74431 - 3.223 n^2) Cos[3.14159 n] +
5.79563 n Sin[3.14159 n])/(n (6.29339 + 1. n^2)) +
1/(n (5.35495 + 1. n^2)) ((4.75667 + 3.27006 n^2) Cos[
3.14159 n] + (-4.75667 - 0.889934 n^2) Cos[6.28319 n] +
5.51162 n Sin[3.14159 n] - 0.00383716 n Sin[6.28319 n])))*
Sin[10 (n) π (x)], {n, 1, 100}]


The graph of this fit (100 partial sums):

Finally I used Mathematica to evaluate the PDE here is my code:

κ = 2.61*10^-3;
α = 1.12*10^-4;
NDSolve[{
Derivative[0, 1][u][x, t] == α*Derivative[2, 0][u][x, t] - κ*(u[x, t] - 17),
Derivative[1, 0][u][0, t] == 0,
Derivative[1, 0][u][0.2, t] == 0,
u[x, 0] == (0.7607201329841713/2 +
Sum[(1/10 ((0.912752572216549 Sin[3.141592653589793 n])/
n + (-0.002189478133985855 +
5.795634757369874 Cos[3.141592653589793 n] +
2.310248484743582 n Sin[3.141592653589793 n])/(
6.293389749792516 +
1. n^2) + (-0.888276236822269 Sin[
3.141592653589793 n] +
0.888276236822269 Sin[6.283185307179586 n])/n +
1/(5.3549458751438115 +
1. n^2) (5.51162427937814 Cos[3.141592653589793 n] -
0.003837157071497423 Cos[6.283185307179586 n] -
2.3817822905923016 n Sin[3.141592653589793 n] +
0.0016581813809965903 n Sin[6.283185307179586 n])))*
Cos[10 (n) π (x)] + (1/
10 ((5.74431 +
0.913625 n^2 + (-5.74431 - 3.223 n^2) Cos[3.14159 n] +
5.79563 n Sin[3.14159 n])/(n (6.29339 + 1. n^2)) +

1/(n (5.35495 + 1. n^2)) ((4.75667 + 3.27006 n^2) Cos[
3.14159 n] + (-4.75667 - 0.889934 n^2) Cos[
6.28319 n] + 5.51162 n Sin[3.14159 n] -
0.00383716 n Sin[6.28319 n])))*
Sin[10 (n) π (x)], {n, 1, 100}])},
u, {t, 0, 10}, {x, 0, 0.2}]


However the plot of this solution using Mathematica, it gives me a very funny graph:

This is completely wrong as the bar is supposed to loose heat not gain heat. Any ideas as to what is going on? And the PDE comes with this error:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent. >>

Here is the plot of the actual data, the solution to the PDE should look something like this:

Any help will much appreciated.

• Your definition of u[x, 0] is definitely inconsistent with your other equations. Given your expression for u[x, 0], I get 1.67605 for Derivative[1, 0][u][0., 0] and 0.198103 for Derivative[1, 0][u][0.2, 0]. So the error message you are getting seems cogent, and your first task is to specify a consistent system of equations. – m_goldberg Jul 5 '15 at 19:35
• Inconsistent boundary conditions certainly is undesirable but is unlikely to be causing u to grow with time, because this problem is localized to small t. Instead, - κ*(u[x, t] - 17) appears to be the problem, acting as a source of heat as long as u < 17. If this term is meant to represent lateral heat loss, perhaps it should be part of the boundary conditions in x. Even then, it still may be a source of heat until for u < 17. – bbgodfrey Jul 5 '15 at 22:27
• As for the boundary conditions, they are zero because the ends of the bar are insulated. Therefore the spatial derivative (heat flux) at the 2 ends are zero. But how do I make a function that would be consistent with the initial temperature distribution? – daniel17027 Jul 6 '15 at 1:00