Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

RE: Re: PDE heat equation (inconsisten problem)

• To: mathgroup at smc.vnet.net
• Subject: [mg93340] RE: [mg93233] Re: [mg93177] PDE heat equation (inconsisten problem)
• From: "Ingolf Dahl" <ingolf.dahl at telia.com>
• Date: Tue, 4 Nov 2008 06:18:41 -0500 (EST)
• Organization: University of Gothenburg
• References: <200810300701.CAA00581@smc.vnet.net> <200811011004.FAA15491@smc.vnet.net>
• Reply-to: <ingolf.dahl at physics.gu.se>

Hi Matteo,
Your question is closely related to my question in
http://forums.wolfram.com/mathgroup/archive/2008/Aug/msg00532.html 
and the following discussion. I have not yet said my last word in that
discussion.

I think that your boundary and initial conditions in x=0, t=0 are correct and
consistent T[0,0]=40, contrary to what Tony Harker is stating. The initial
condition is discontinuous, and I think that that is OK mathematically, even
if Mathematica has trouble to represent this function as an interpolation
function.

The problem is instead located at x=side (where side=1) and t=0. If you
evaluate the boundary condition 

bc2 = {Derivative[1, 0][T][side, t] == -h (T[side, t] - to)}

at t=0, you will find that the initial condition implicates that
Derivative[1, 0][T][side, t] = 0, while -h (T[side, t] - to)= 20. We thus
have an apparent inconsistency.
If we change your initial condition to 

ic = {T[x, 0] == If[x > 0, 20., 40]};

there will be no inconsistency at x=side, t=0. The following warning 

"NDSolve::ibcinc: Warning: Boundary and initial conditions are inconsistent.
>>" 

vanishes, and we obtain an credible solution. If we compare this solution
with your initial solution, we see that your initial solution tends a
steady-state solution where T[side, t] tends to 20, while the new solution
here tends to 30, which seems more correct. This tells us that Mathematica
is not able to cope with this kind of "inconsistency" in the correct way.
What really happens in a physical system is that the system in very short
time will adapt to the boundary condition bc2. At the time t=0 we should not
use the initial condition ic to calculate the first order derivative
Derivative[1, 0][T][side, t]. The derivative Derivative[1, 0][T][x,t]
instead is discontinuous at the border, and just at the border at time t=0
it should take the value determined by bc2. I have not been able to force
NDSolve and Mathematica to perform this trick, and I consider this as a kind
of bug. Other people might have another opinion about that. 

However, we might modify your system a little, to allow a slower adoption to
the boundary condition bc2. We put a thin metal layer on the end of your 
silicon bar, with a finite non-zero heat capacitance. The layer should have 
enough heat conductance to be able to keep at one uniform temperature
T[side, t]. 
By letting the heat capacitance of this layer approach zero, we should 
obtain your initial problem again.
We might then modify your equations in the following way

kappa = 1.;
ro = 1.;
c = 1.;
k = (ro*c)/kappa;
h = 1.;
side = 1.;(*silicon bar length*)
tmax = 2;
to = 20.;(*room temperature*)
alfa = 0.01; (* corresponds to heat capacity of the metal layer *)
eq1 = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[x, t]\)\) == 
   Which[x == 0, 0., 
    x == side, -(Derivative[1, 0][T][x, t] + h (T[x, t] - to))/alfa, 
    True, Derivative[2, 0][T][x, t]/k]}; 
ic = {T[x, 0] == If[x > 0, 0., 40.]}; 
bc1 = {T[0, t] == 40.};
sol1 = NDSolve[{Join[eq1, ic, bc1]}, 
  T[x, t], {x, 0, side}, {t, 0, tmax}];

Note how I have move the boundary condition bc2 into the differential
equation.
Mathematic grumbles a warning "NDSolve::bcart: Warning: An insufficient
number of
boundary conditions have been specified for the direction of independent
variable y. Artificial boundary effects may be present in the solution",
but delivers a solution that looks appropriate. I have also tested to put 
alfa=0.00000000001 without any problems except the same warning, and
conclude 
that your initial problem must have been correctly specified, but does not 
work due a deficiency of NDSolve.

Best regards

Ingolf Dahl
Sweden
ingolf.dahl at telia.com

> -----Original Message-----
> From: Tony Harker [mailto:a.harker at ucl.ac.uk] 
> Sent: den 1 november 2008 11:04
> To: mathgroup at smc.vnet.net
> Subject: [mg93233] Re: [mg93177] PDE heat equation 
> (inconsisten problem)
> 
> 
> What Mathematica is telling you is quite right. Your initial 
> conditions implies T[0,0]=0, and your boundary condition bc1 
> implies T[0,0]=40, so they are incompatible. This is a fairly 
> standard situation, and physicists and engineers annoy 
> mathematicians with their rather cavalier attitude to it, but 
> your answer is probably fine.
> 
>   Tony Harker
> 
> Dr A.H. Harker
> Department of Physics and Astronomy
> University College London
> Gower Street
> London
> WC1E 6BT
> 
> Tel: (44)(0) 2076793404
> E:    a.harker at ucl.ac.uk
> 
>  EDUCATION, n. That which discloses to the wise and disguises 
> from the foolish their lack of understanding. (Ambrose 
> Bierce, The Devil's Dictionary, 1911)
> 
> 
> ]-> -----Original Message-----
> ]-> From: Matteo Calabrese [mailto:calabrex87 at hotmail.it] ]-> 
> Sent: 30 October 2008 07:01 ]-> To: mathgroup at smc.vnet.net 
> ]-> Subject: [mg93177] PDE heat equation (inconsisten 
> problem) ]-> ]-> Dear Mathematica Friends, ]-> ]-> I'm trying 
> to solve this simple problem: I've got a silicon ]-> bar in 
> 1D resolving fourier equation of heat. With these ]-> 
> Boundary Condition, mathematica gives me this kind of ]-> 
> error. however, solution seems consistent with the problem.
> ]-> Could anybody explain ad resolve it??
> ]->
> ]-> thank you in advance,
> ]->
> ]-> Matteo Calabrese
> ]->
> ]-> University of Physics
> ]->
> ]-> Turin, Italy
> ]->
> ]->
> ]-> (*Mathematica Vs 6.0*)
> ]->
> ]->
> ]->
> ]-> kappa=1.;
> ]->
> ]-> ro=1.;
> ]->
> ]-> c=1.;
> ]->
> ]-> k=(ro*c)/kappa;
> ]->
> ]-> h=1.;
> ]->
> ]-> side=1.; (*silicon bar length *)
> ]->
> ]-> tmax=1;
> ]->
> ]-> to.; (* room temperature*)
> ]->
> ]->
> ]->
> ]-> eq={Tx,x[x,t]==k*Tt[x,t]}; (*Fourier's heat equation 1D*) 
> ]-> ]-> (*Initial condition*) ]-> ]->  ic={T[x,0]== 
> If[x>0,0,40]}; ]-> ]-> (*Boundary Condition*) ]-> ]->  
> bc1={T[0,t]== 40.}; (*dirichlet condition*) ]-> ]->  
> bc2={Derivative[1,0][T][side,t]== -h (T[side,t]-to)} ]-> 
> (*Newton-Robin condition*) ]-> ]-> ]-> ]->  
> sol=NDSolve[{Join[eq,ic,bc1,bc2]},T[x,t],{x,0,side},{t,0,tmax}]
> ]->
> ]->  NDSolve::ibcinc: Warning: Boundary and initial 
> conditions ]-> are inconsistent. =E2=80=A1 ]-> ]->  = 
> {{T[x,t]=C2=A2=C3=A7InterpolatingFunction[{{0.,1.},{0.,1.}},<>][x,t]}}
> ]->
> ]->
> ]->
> ]-> (*If you plot the Interpolating function, it seems a good 
> solution*) ]-> ]-> ]-> ]-> ]-> ]-> ________ ]->
> 
> 
> 

• References:
  • Re: PDE heat equation (inconsisten problem)
    • From: "Tony Harker" <a.harker@ucl.ac.uk>