PDE problem with Wentzell-Robin boundary condition: wrong solution without any warning
- To: mathgroup at smc.vnet.net
- Subject: [mg93505] PDE problem with Wentzell-Robin boundary condition: wrong solution without any warning
- From: "Ingolf Dahl" <f9aid at chalmers.se>
- Date: Wed, 12 Nov 2008 06:46:35 -0500 (EST)
- Organization: University of Gothenburg
- Reply-to: <ingolf.dahl at physics.gu.se>
Dear Mathematica Friends,
October 30 this year, Matteo Calabrese posted a question, (PDE heat equation
(inconsisten problem)), which I already have answered. Since I am
collaborating with Giovanni Barbero on a related liquid crystal problem, I
have looked a bit further. In the problem Matteo Calabrese posted, NDSolve
claimed inconsistency between the initial and the boundary conditions. I got
the idea to present to NDSolve nice consistent initial and boundary
conditions, starting at time t=-1. Then I let the system evolve, and check
what will happen if I change boundary conditions at time t=0. Will NDSolve
still claim inconsistency? At time t=0 there should not be any difference
between the case that I have let the system evolve since t=-1, or the case
that I specify a new initial condition to be the result of the evolution
since t=-1. Thus this could shed some light upon the initial-boundary
inconsistency.
My approach to this problem is that of a mathematical experimentalist. Maybe
some mathematical theorist in MathGroup could provide more information about
this problem. I hope so. My feeling for the problem is that it maybe was
solved already 1850 or at least 1920.
Imagine a rod, one end kept constant at temperature zero, the other end in
contact with a surrounding media.
First I tried with a smooth transition around t=0, in the following way
te[t_, p_] := 10. + 10*(2/Pi)*ArcTan[p*Tan[t*Pi/2]];
(* te is surrounding temperature, changing from zero at t==-1 to 20. at
t==1, and p is a smoothness parameter, giving sharp transitions for high
values of p. *)
p =.; eq = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(T[x, t]\)\) == \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[x,t]\)\)};
(*Fourier's heat equation 1D*)
(*Initial condition*)
ic = {T[x, -1] == 0.};
(*Boundary Condition*)
bc1 = {T[0, t] == 0.};
(*dirichlet condition*)
bc2 = {Derivative[1, 0][T][1., t] == - (T[1., t] - te[t, p])};
(* Wentzell-Robin condition, is this the correct name? *)
p=10000; (* For this p value, the temperature transition is quite sharp *)
sol1=NDSolve[{Join[eq,ic,bc1,bc2]},T[x,t],{x,0,1},{t,-1,1}];
(* This evaluates without warning. We might plot the solution, and will then
see a plausible solution *)
Plot3D[T[x, t] /. sol1, {x, 0, 1}, {t, -1, 1}]
(* We can see that the x-derivative of T[x,t] is almost discontinuous at
x==1, t==0 *)
Needs["NumericalCalculus`"]
Plot3D[ND[T[x, t] /. sol1, x, x0], {x0, 0, 1}, {t, -1, 1}]
Now we make a sharp transition in surrounding temperature instead:
te[t_] := 20*UnitStep[t];
eq = {\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(T[x, t]\)\) == \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[x,t]\)\)};
(*Fourier's heat equation 1D*)
(*Initial condition*)
ic = {T[x, -1] == 0.};
(*Boundary Condition*)
bc1 = {T[0, t] == 0.};
(*dirichlet condition*)
bc2 = {Derivative[1, 0][T][1., t] == - (T[1., t] - te[t])};
sol2 = NDSolve[{Join[eq, ic, bc1, bc2]},
T[x, t], {x, 0, 1}, {t, -1, 1}];
Plot3D[T[x, t] /. sol2, {x, 0, 1}, {t, -1, 1}]
The evaluation is without any warning, but the plot shows a function, which
is completely flat and evidently very wrong!
It is reasonable to assume that the two solutions should be similar. I think
that Mathematica here completely misses the discontinuity of the
x-derivative of T[x,t] at x==1, t==0 , and that the problems Matteo
Calabrese had was related to this. Any comments?
Ingolf Dahl
Sweden