       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

```Dear Mathematica Friends,

October 30 this year, Matteo Calabrese posted a question, (PDE heat equation
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
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

Ingolf Dahl

Sweden

```

• Prev by Date: Re: Storing and Loading Definitions / Emulating Associative Arrays
• Next by Date: Re: Re: How do little quickest ?
• Previous by thread: Presentations Solutions
• Next by thread: Visulizing the content of a variable