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

**Re: Storing and Loading Definitions / Emulating Associative Arrays**

**Re: Re: How do little quickest ?**

**Presentations Solutions**

**Visulizing the content of a variable**