PDE heat diffusion to sphere
- To: mathgroup at smc.vnet.net
- Subject: [mg7735] PDE heat diffusion to sphere
- From: "w.meeussen" <meeussen.vdmcc at vandemoortele.be>
- Date: Thu, 3 Jul 1997 01:28:43 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
before I go off on vacation, a query for suggentions on how to clean up my
Mma 3.0 syntax for a problem that "just by the skin of the teeth" works out
right.
I'v got the solution,
I still need the most economical & elegant way to re-produce & generate it.
****************************************************************************
*************
Problem:
Unsteady state heat conduction from a infinite bulk towards a sphere.
Boundary conditions:
Temperature T at infty is 0; T on sphere (radius b) is 1 at all times (t):
****************************************************************************
*************
the un-elegant way :
<<Calculus`VectorAnalysis`
SetCoordinates[Spherical[r, th, fi]]
eqn1=dif Laplacian[f[r,t]] == D[f[r,t],t]
<<Calculus`LaplaceTransform`
LaplaceTransform[f[r_,t_],t_,s_]:=U[s][r]
eqn2=LaplaceTransform[eqn1,t,s]
bound={U[s][d]==1/s,U[s][\[Infinity]]==0}
(* boundary conditions are a real headache *)
it=Join[{eqn2/.f[r,0]->0},bound]
DSolve[it,U[s][r],r]
(* produces warnings & fat output including constant
C[2] *)
{{U[s][r] -> (2*d*E^(d*Sqrt[s/dif]) -
(E^(2*d*Sqrt[s/dif]) - E^(2*r*Sqrt[s/dif]))*Sqrt[s/dif]*C[2])/
(2*E^(r*Sqrt[s/dif])*r*s)}}
(* hoping for the best : put C[2] to zero and touch
wood *)
%/.C[2]->0
{{U[s][r] -> (d*dif*E^(d*Sqrt[s/dif] - r*Sqrt[s/dif] - Log[dif*r]))/s}}
(* going for it ... *)
InverseLaplaceTransform[%,s,t,Assumptions->{Re[(d-r)^2/dif]>0,t>0}]
(* this produces something good, packed in FAT If[
.. ]'s *)
(* extracting the good part : *)
d*dif*
(dif^((3*-1)/2)*(d - r)*(d - r)^((2*-1)/2)*(-1 + Erf[(d - r)^(2/2)/
(Sqrt[dif]*2*Sqrt[t])]))/(1/Sqrt[dif])/r
(* massaging produces a good solution: *)
f[r_,t_]:=1-((d*Erfc[-(d - r)/(2*Sqrt[dif]*Sqrt[t])])/r)
check it by letting the PDE "eqn1" loose on it :
eqn1//Simplify gives "True".
The limit for t->infty looks like it should:
Limit[f[r,t],t->\[Infinity]] gives: 1-d/r
conclusion: f[r,t] is a viable solution to the physics problem.
****************************************************************************
********
Mma 3.0 problems: the InverseLaplaceTransform eats up more than my 32 Mbyte
RAM if I don't break up the problem by calculating & saving intermediates,
then restarting Mma from there.
Although I use Assumptions->{Re[(d-r)^2/dif]>0,t>0}
the output still contains an If[..] like:
If[t > 0 && Re[(d - r)^2/dif] > 0 &&
Re[-(d^2/(4*dif)) + (d*r)/(2*dif) - r^2/(4*dif)] < 0, ...
which seems un-logical (?)
How should I define boundary conditions for something "un-physical" like U[x][r]
in order to avoid getting a solution containing "C[2]" ?
Anyone with an idea what the InverseLaplaceTransform of the "entire"
solution would look like (not "arbitrarily" setting C[2]->0)?
Yes, I realise that it's risky medling like this with insufficient math
background,
I promis a will by a book one of these days.
If it weren't for the checks using the PDE "eqn1" and
Limit[f[r,t],t->infty], I wouldn't trust it. But how should it be done right?
take your time for answering this,
I'm gone fishing for three weeks,
wouter.
Dr. Wouter L. J. MEEUSSEN
eu000949 at pophost.eunet.be
w.meeussen.vdmcc at vandemoortele.be