Re: NDSolve
- To: mathgroup at smc.vnet.net
- Subject: [mg91400] Re: NDSolve
- From: Jean-Marc Gulliet <jeanmarc.gulliet at gmail.com>
- Date: Wed, 20 Aug 2008 04:04:25 -0400 (EDT)
- Organization: The Open University, Milton Keynes, UK
- References: <g8e9sn$48n$1@smc.vnet.net>
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] )
ldrp[x,a,b]:=kg-(kg-kp)*heat[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.
In[11]:=
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
heat[x,5.*10^-6,1.*10^-6]+ldrp[x,5.*10^-6,1.*10^-6]
(u^(0,2))[t,x],u[0,x]==0,True,u[t,201/200000]==0},u,
{t,0,tmax},{x,-(201/200000),201/200000},
MaxSteps->25000,MaxStepFraction->0.0002,PrecisionGoal->2]
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.
Regards,
-- Jean-Marc