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