PDE problem with Wentzell-Robin boundary condition: wrong solution without any warning

Discussion in 'Mathematica' started by Ingolf Dahl, Nov 12, 2008.

  1. Ingolf Dahl

    Ingolf Dahl Guest

    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
     
    Ingolf Dahl, Nov 12, 2008
    #1
    1. Advertisements

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.