Re: diffusion PDE
- To: mathgroup at smc.vnet.net
- Subject: [mg131021] Re: diffusion PDE
- From: sabri <sabrinacasanova at gmail.com>
- Date: Tue, 4 Jun 2013 01:59:43 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- Delivered-to: l-mathgroup@wolfram.com
- Delivered-to: mathgroup-outx@smc.vnet.net
- Delivered-to: mathgroup-newsendx@smc.vnet.net
- References: <koch1r$t5b$1@smc.vnet.net> <koeh9t$2rv$1@smc.vnet.net>
Hello Roland,
Thanks for your reply. I thought to have imposed the right condition, i.e. startcond = 0 everywhere except in at the boundary for radius = radiusmax. But I must have made a syntax error. Yes, you are right. I would expect the start condition to be a steady state solution. But the solution should come to the steady state through intermediate states.
When I put as startcond = 0 (see notebook below), then I still get an error (see below) and I do not get the steady state solution I would expect. Are you sure that I can neglect the error message ?
Best
Sabrina
------------------------
chi = 1.;
d0 = 3. 10^27;
n = 300.;
norm = 4. 3.14 1.8/(3 10^10);
B = 100. Sqrt[n/10000.];
Dif = chi*d0/Sqrt[B/3];
tau0 = 2.*10^5/(n/300.)*3.10^7;
gam = 0.5;
radiusmin = 10^(-5);
radiusrad = 20.;
radiusmax = 40.;
timemin = 0.;
timemax = 1000 (40.*3.18 10^18)^2/(6.*Dif)
energymin = 1.;
energymax = 10^8;
energycut = 10^6;
pde = D[y[time, radius, energy], time] -
Dif*energy^gam*D[y[time, radius, energy], radius, radius] -
Dif*energy^gam/radius*2*D[y[time, radius, energy], radius] +
y[time, radius, energy]/tau0 +
energy D[y[time, radius, energy], energy]/tau0 == 0
startcond = y[timemin, radius, energy] == 0
rmaxcond =
y[time, radiusmax, energy] == norm/energy^2.7*Exp[-energy/energycut]
rmincond = Derivative[0, 1, 0][y][time, radiusmin, energy] == 0
energymaxcond = y[time, radius, energymax] == 0
energymincond =
y[time, radius, energymin] == startcond[[2]] /. energy -> energymin
s = NDSolve[{pde, startcond, rmincond, rmaxcond, energymaxcond},
y, {radius, radiusmin, radiusmax}, {energy, energymin,
energymax}, {time, timemin, timemax},
MaxStepFraction -> {1/1000, 1/40, 1/40}]
Plot3D[y[time, radiusmax, energy] /. s, {time, timemin,
timemax}, {energy, energymin, energymax/10}, PlotRange -> All,
PlotPoints -> 200, AxesLabel -> {"time", "energy", "y"}]
LogLogPlot[
Evaluate[y[timemax, radiusmax, energy] /. s], {energy, energymin,
energymax}]
NDSolve::eerr: Warning: scaled local spatial error estimate of 1.2872504795358698`*^16 at time = 2.1598391239374475`*^15 in the direction of independent variable radius is much greater than the prescribed error tolerance. Grid spacing with 41 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options. >>
---------------------------------------------------
Sunday, June 2, 2013 6:22:21 AM UTC+2, Roland Franzius wrote:
> Am 01.06.2013 12:05, schrieb sabrinacasanova at gmail.com:
>
> > Hello,
>
> >
>
> >
>
> > In the notebook below you can find PDE which represents a particle diffusion equation
>
> > with energy dependence and not only time/spatial dependences. The
>
> > particles should diffuse in a energy dependent way inside, lets say, a
>
> > sphere. There is a spherical symmetry, so only the radius of this
>
> > sphere is important. I have an initial condition and a boundary
>
> > condition on the spatial border of the sphere. My question is whether Mathematica
>
> > is able to solve such PDEs at all. If yes, do you know of somebody
>
> > who has already done it ?
>
> > Best
>
> >
>
> > Sabrina
>
> >
>
> > --------------------------------
>
> >
>
> > chi = 1.
>
> > d0 = 3. Power[10., 27]
>
> > n = 300.
>
> >
>
> > norm = 4. 3.14 1.8/( 3 Power[10., 10.])
>
> > Evaluate [norm]
>
> > B = 100. Power[n/10000., 0.5]
>
> > Dif = chi*d0*Power[B/3., -0.5]
>
> > tau0 = 2.*Power[10., 5]*Power[(n/300.), -1] 3. Power[10., 7]
>
> > gam = 0.5
>
> >
>
> > x2min = 0.00001
>
> > x2rad = 20.
>
> > x2max = 40.
>
> > x1min = 0.
>
> > x1max = Power[(40.*3.18 Power[10., 18.]), 2]/(6.*Dif)
>
> > x3min = 1.
>
> > x3max = 100000000000.
>
> > x3cut = 1000000
>
> >
>
> > pde = D[y[x1, x2, x3], x1] -
>
> > Dif Power[x3, gam] D[y[x1, x2, x3], {x2, 2}] -
>
> > Dif Power[x3, gam] Power[x2, -1] 2 D[y[x1, x2, x3], x2] +
>
> > y[x1, x2, x3]/tau0 + x3 D[y[x1, x2, x3], x3]/tau0 == 0
>
> >
>
> >
>
> > s = NDSolve[{pde,
>
> > y[x1min, x2, x3] ==
>
> > norm Power[x3, -2.7] Exp[-x3/x3cut] UnitStep[x2 - x2max] UnitStep[
>
> > x3max - x3],
>
> > y[x1, x2max, x3] ==
>
> > norm Power[x3, -2.7] Exp[-x3/x3cut] UnitStep[x3max - x3],
>
> > Derivative[0, 1, 0][y][x1, x2min, x3] == 0,
>
> > y[x1, x2, x3max] == 0 },
>
> > y, {x1, x1min, x1max}, {x2, x2min, x2max}, {x3, x3min, x3max}]
>
> >
>
> > Plot3D[Evaluate[y[x1, x2max, x3] /. s], {x1, x1min, x1max}, {x3,
>
> > x3min, x3max}, PlotRange -> All]
>
> >
>
>
>
>
>
>
>
> I've written the system a bit more readable.
>
>
>
> There seems to be no problem exept a lacking boundary condition at
>
> energymin and the fact, that the solution seems to be independent of
>
> time. The start condition seems to be already a solution, make a test of=
>
> it applying the operator pde on the start condition.
>
>
>
>
>
> chi = 1.;
>
> d0 = 3. 10^27;
>
> n = 300.;
>
> norm = 4. 3.14 1.8/(3 10^10);
>
> B = 100. Sqrt[n/10000.];
>
> Dif = chi*d0/Sqrt[B/3];
>
> tau0 = 2.*10^5/(n/300.)* 3.10^7;
>
> gam = 0.5;
>
> radiusmin = 0.00001;
>
> radiusrad = 20.;
>
> radiusmax = 40.;
>
> timemin = 0.;
>
> timemax = 1000 (40.*3.18 10^18)^2/(6.*Dif);
>
> energymin = 1.;
>
> energymax = 100000000000.;
>
> energycut = 1000000;
>
> pde = D[y[time, radius, energy], time] -
>
> Dif *energy^gam *D[y[time, radius, energy], radius, radius] -
>
> Dif* energy^gam/radius * 2* D[y[time, radius, energy], radius] +
>
> y[time, radius, energy]/tau0 +
>
> energy D[y[time, radius, energy], energy]/tau0 == 0
>
>
>
> startcond =
>
> y[timemin, radius, energy] == norm/ energy^2.7 Exp[-energy/energycu=
t]
>
>
>
> rmaxcond =
>
> y[time, radiusmax, energy] == norm/energy^2.7* Exp[-energy/energycu=
t]
>
> rmincond = Derivative[0, 1, 0][y][time, radiusmin, energy] == 0
>
> energymaxcond = y[time, radius, energymax] == 0
>
> energymincond =
>
> y[time, radius, energymin] == startcond[[2]] /. energy -> energymin
>
>
>
>
>
>
>
> s = NDSolve[{pde, startcond, rmincond, rmaxcond, energymincond,
>
> energymaxcond},
>
> y, {radius, radiusmin, radiusmax}, {energy, energymin,
>
> energymax}, {time, timemin, timemax},
>
> MaxStepFraction -> {1/1000, 1/40, 1/40}]
>
>
>
>
>
> Plot3D[y[time, radiusmax, energy] /. s, {time, timemin,
>
> timemax}, {energy, energymin, energymax/10}, PlotRange -> All,
>
> PlotPoints -> 200, AxesLabel -> {"time", "energy", "y"}]
>
>
>
> --
>
>
>
> Roland Franzius