MathGroup Archive 2008

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

Search the Archive

Re: NDSolve

  • To: mathgroup at
  • Subject: [mg91400] Re: NDSolve
  • From: Jean-Marc Gulliet <jeanmarc.gulliet at>
  • Date: Wed, 20 Aug 2008 04:04:25 -0400 (EDT)
  • Organization: The Open University, Milton Keynes, UK
  • References: <g8e9sn$48n$>

Uli Wuerfel wrote:

> I was trying to solve the heat equation with NDSolve(one spatial dimension, x, and time t). The simulated device (centered around x=0) consists of three regions: Two rather big volumina (with thickness each being dg) to the left and right of a thin layer (with thickness dp). A heat source is only present in the thin layer. The outer left boundary of the whole device is open, i.e. the spatial derivative of the temperature is zero there. At the right boundary, there is perfect cooling, i.e. the temperature stays constant. 
> As the heat is generated only in the thin layer, the change is quite abrupt. What I would like to know is if there is a way to tell Mathematica to place a very fine grid there because the result I get is not correct.
> Here is the code:
> dp=1e-5;
> dg=1e-3;
> eta=500/dp;
> sol=u/.First[
> NDSolve[D[u[t, x],t]==ldrp[x,0.5*dp,0.1*dp]*D[u[t,x],(x,2)]+eta*heat[x,0.5*dp,0.1*dp],
> u[0,x]==0,D[u[t,-dg-dp/2],x]==0,u[t,dg+dp/2]==0},u,{t,0,tmax},{x,-dg-dp/2,dg+dp/2},MaxSteps->25000,MaxStepFraction->0.0002,PrecisionGoal->2]]
> with
> heat[x,a,b]:=-1+1/(Exp[(x - a)/b]+1)+1/(Exp[-(x+a)/b]
> being the (steep) function describing the extension of the heat source and
> ldrp[x,a,b]:=kg-(kg-kp)*heat[x,a,b]
> describing the spatial characteristics of the heat conductivity divided by the product of mass density and heat capacity.
> The result depends unfortunately very much on the thickness of the thin layer, although the TOTAL heat generated there is always the same (as eta=500/dp).
> Also, the steady-state temperature distribution in the right volume (with the cooling at the outer right boundary) where there is no heat generation should be
> just a function of the total generated heat and the heat conductivity of that volume. So the temperature difference
> should (in the steady-state) always be: 
> u[x=dp/2]-u[x=dg+dp/2]=500/lambda_g, 
> with lambda_g being the heat conductivity of that volume.
> But the results produced by Mathematica are different.
> Do you know what I am doing wrong or is it just that it
> is not possible to handle such abrupt changes?
> Thank you for any advice.

Mathematica has few syntactic rules and convention; however, these rules 
and conventions must be followed rigorously and they may differ 
significantly from other programming languages and CAS.

For instance, the following input is almost certainly not what you want:

   In[1]:= dp = 1e-5
           dg = 1e-3
           eta = 500/dp

   Out[1]= -5+e
   Out[2]= -3+e
   Out[3]= 500/(-5+e)

I believe that what you meant is

   In[4]:= dp = 10^-5
           dg = 10^-3
           eta = 500/dp

   Out[4]= 1/100000
   Out[5]= 1/1000
   Out[6]= 50000000

Also, parentheses and curly brackets are not interchangeable:

   In[12]:= D[u[t,x],(x,2)]

     During evaluation of In[12]:= Syntax::sntxf: "(" cannot be
      followed by "x,2)".

     During evaluation of In[12]:= Syntax::tsntxi: "x,2" is
      incomplete; more input is needed.

     During evaluation of In[12]:= Syntax::sntxi: Incomplete
      expression; more input is needed.

The correct syntax is:

   In[7]:= D[u[t,x],{x,2}]

   Out[7]= (u^(0,2))[t,x]

The following partial derivative always evaluates to zero since it is 
independent of the variable x:

   In[8]:= D[u[t,-dg-dp/2],x]

   Out[8]= 0

   In[9]:= heat[x,a,b]:=-1+1/(Exp[(x-a)/b]+1)+1/(Exp[-(x+a)/b] )

Note that NDSolve failed because you have not provided a value for tmax 
and the condition D[u[t, -dg - dp/2], x] == 0 is meaningless for it is 
always true.

    NDSolve[{D[u[t, x], t] ==
      ldrp[x, 0.5*dp, 0.1*dp]*D[u[t, x], {x, 2}] +
       eta*heat[x, 0.5*dp, 0.1*dp], u[0, x] == 0,
     D[u[t, -dg - dp/2], x] == 0, u[t, dg + dp/2] == 0}, u, {t, 0,
     tmax}, {x, -dg - dp/2, dg + dp/2}, MaxSteps -> 25000,
    MaxStepFraction -> 0.0002, PrecisionGoal -> 2]

     During evaluation of In[11]:= NDSolve::deqn: Equation or
      list of equations expected instead of True in the
      first argument {<<1>>}.

   Out[11]= NDSolve[{(u^(1,0))[t,x]==50000000

So, you should post some code on MathGroup (this list/newsgroup) that 
one can easily cut and past into Mathematica and that reproduces the 
problem you have.

-- Jean-Marc

  • Prev by Date: Re: Need a Faster Solution to Number Theory Problem
  • Next by Date: Re: Conflicting Switch Statements
  • Previous by thread: NDSolve
  • Next by thread: Fourier-Motzkin elimination