meditationatae

Just another WordPress.com site

numerical methods, heat flow

For a heated wire of length 2pi which
is insulated against heat flow at
both ends, we can model the wire
and its temperature distribution
by a function on the segment
[0, 2pi].

We can suppose that a point heat source
is at some point `a’ between the two
extremities.  In the example below,
a = 1.  In physical terms, at time
t = 0, all of the wire is at the
temperature of absolute zero, except
for a Dirac delta pseudo-function
representing a source of heat of
infinite temperature and non-zero
heat content at `a’.

As time goes by, heat diffuses through
the wire while the average temperature
remains connstant over the whole length
of 2pi because heat can neither leave
nor escape the wire, the segment of
length 2pi.

As time goes by, the heat distribution
is more and more uniform.  The point
`a’ is closer to the end at 0 than to
the end at 2pi, since we set a = 1.
Therefore, we can expect point 0 to be
at a higher temperature than point 2pi.

One can attempt to solve the heat equation
here by separation of variables.  As is
often the case, this leads to eigenvalue
problems for differential operators.
In this case, the differential operator
is the Laplacian operator with Neumann
boundary conditions.

The Neumann boundary conditions are
equivalent to requiring of a
potential eigenfunction u that
u'(0) = 0 and u'(2pi) = 0.
In other words, the gradient of u
is zero at both heat-insulated ends
of  the wire.  

For an animation, see for example:
“Solution of a 1D heat equation PDE”
at Wikipedia:
http://en.wikipedia.org/wiki/Heat_equation

In sci.math in the last few days,
I sketched an attempt at a numerical
solution based on “mirror images”
of the standard heat kernel as
one “unfolds” the domain accross
the boundary points (or edges of
a triangle, square) and then repeats
the “unfolding” over and over as
with billiards problems for
tables that have a triangular or
rectangular boundary.

Cf.: Discussion thread including this post:
http://mathforum.org/kb/message.jspa?messageID=7930875

In Pari/GP, I gave function definitions represented by:

func =
  (Q)->sum(N=0,maxi,ker(Q,a+4*Pi*N))+sum(N=0,maxi,ker(Q,(2*Pi-a)+2*Pi+4*Pi*N))+sum(N=1,maxi,ker(Q,a-4*Pi*N))+sum(N=1,maxi,ker(Q,(2*Pi-a)+2*Pi-4*Pi*N))

ker =
  (X,W)->exp(-(X-W)*(X-W)/T)/sqrt(T)

shaq =
  (Q)->(1/2)+cos((Q/50)*Pi)/2

with constants:
a = 1,    T = 400 and  maxi = 150 .

T is a time-value of 400 units after
time = 0.  By that time, the temperature
distribution is very nearly uniform.

The point source is located at a = 1.

The parameter maxi reflects the fact that
about 4*maxi = 600 intervals of length 2pi
with their “image sources” are counted.
Images appear because the ends of the wire
at 0 and 2pi are insulated, and heat that’s
reflected is similar to heat that goes over
to the image through the boundary (as with
billiards problems).

‘ker’ is the heat kernel function at time T.
I imagine the source at formal parameter W
and the evaluation point being the formal
parameter X.

‘func’ is the implementation of the image method
for this problem.  There is a sum of about
500 Gaussian terms.

‘shaq’ is a normalization of the
x |-> cos(x/2)  eigenfunction
with x in [0, 2pi] dilated to
[0, 50] and the codomain [-1, 1]
mapped affinely to [0, 1].

The Gaussian heat kernel for Euclidian spaces
is explained here:
http://en.wikipedia.org/wiki/Heat_kernel

Then, I find that at T=400 seconds, the
normalized numerical solution given by
‘func’ agrees with the exact eigenvalue
represented by ‘shaq’ to about
26 decimal digits.

Some PARI/gp output is given below:

? for(m=0,50,print(m,”  “,(func((m/50)*2*Pi)-func(2*Pi))/(func(0)-func(2*Pi))))
0   1.000000000000000000000000000
1   0.999013364214135780976168403
2   0.996057350657238915524896521
3   0.991143625364344340542820871
4   0.984291580564315559745084187
5   0.975528258147576786058219666
6   0.964888242944125701830471278
7   0.952413526233009763856834323
8   0.938153340021931793654057951
9   0.922163962751007539274279031
10  0.904508497187473712051146708
11  0.885256621387894615401504818
12  0.864484313710705761573365159
13  0.842273552964344336866141678
14  0.818711994874344855088356405
15  0.793892626146236564584352977
16  0.767913397489498309135654383
17  0.740876837050857637493595751
18  0.712889645782536324431251222
19  0.684062276342338979578473573
20  0.654508497187473712051146708
21  0.624344943582427394121141872
22  0.593690657292862315271275367
23  0.562666616782152122686559379
24  0.531395259764656688038089112
25  0.500000000000000000000000000
26  0.468604740235343311961910887
27  0.437333383217847877313440620
28  0.406309342707137684728724632
29  0.375655056417572605878858127
30  0.345491502812526287948853291
31  0.315937723657661020421526426
32  0.287110354217463675568748777
33  0.259123162949142362506404248
34  0.232086602510501690864345616
35  0.206107373853763435415647022
36  0.181288005125655144911643594
37  0.157726447035655663133858321
38  0.135515686289294238426634840
39  0.114743378612105384598495181
40  0.095491502812526287948853291
41  0.077836037248992460725720967
42  0.061846659978068206345942048
43  0.047586473766990236143165676
44  0.035111757055874298169528722
45  0.024471741852423213941780333
46  0.015708419435684440254915812
47  0.008856374635655659457179128
48  0.003942649342761084475103478
49  0.000986635785864219023831596
50  0.E-27

? for(m=0,50,print(m,”  “,shaq(m)))
0   1.000000000000000000000000000
1   0.999013364214135780976168403
2   0.996057350657238915524896521
3   0.991143625364344340542820871
4   0.984291580564315559745084187
5   0.975528258147576786058219666
6   0.964888242944125701830471278
7   0.952413526233009763856834323
8   0.938153340021931793654057951
9   0.922163962751007539274279031
10  0.904508497187473712051146708
11  0.885256621387894615401504818
12  0.864484313710705761573365159
13  0.842273552964344336866141678
14  0.818711994874344855088356405
15  0.793892626146236564584352977
16  0.767913397489498309135654383
17  0.740876837050857637493595751
18  0.712889645782536324431251222
19  0.684062276342338979578473573
20  0.654508497187473712051146708
21  0.624344943582427394121141873
22  0.593690657292862315271275367
23  0.562666616782152122686559379
24  0.531395259764656688038089112
25  0.500000000000000000000000000
26  0.468604740235343311961910887
27  0.437333383217847877313440620
28  0.406309342707137684728724632
29  0.375655056417572605878858126
30  0.345491502812526287948853291
31  0.315937723657661020421526426
32  0.287110354217463675568748777
33  0.259123162949142362506404248
34  0.232086602510501690864345616
35  0.206107373853763435415647022
36  0.181288005125655144911643594
37  0.157726447035655663133858321
38  0.135515686289294238426634840
39  0.114743378612105384598495181
40  0.095491502812526287948853291
41  0.077836037248992460725720968
42  0.061846659978068206345942048
43  0.047586473766990236143165676
44  0.035111757055874298169528721
45  0.024471741852423213941780333
46  0.015708419435684440254915812
47  0.008856374635655659457179128
48  0.003942649342761084475103478
49  0.000986635785864219023831596
50  0.E-77

end of article.
=========================================
For reference, all that is below.

? T
%126 = 400

? maxi
%127 = 150

? a
%128 = 1

? \u

func =
  (Q)->sum(N=0,maxi,ker(Q,a+4*Pi*N))+sum(N=0,maxi,ker(Q,(2*Pi-a)+2*Pi+4*Pi*N))+sum(N=1,maxi,ker(Q,a-4*Pi*N))+sum(N=1,maxi,ker(Q,(2*Pi-a)+2*Pi-4*Pi*N))

ker =
  (X,W)->exp(-(X-W)*(X-W)/T)/sqrt(T)

shaq =
  (Q)->(1/2)+cos((Q/50)*Pi)/2

Advertisements

Written by meditationatae

December 4, 2012 at 10:36 am

Posted in History

%d bloggers like this: