PDE backward solution problem
- To: mathgroup at smc.vnet.net
- Subject: [mg116930] PDE backward solution problem
- From: Arturas.Acus at tfai.vu.lt
- Date: Sat, 5 Mar 2011 06:03:57 -0500 (EST)
Dear Group,
I am trying to solve ambipolar diffusion PDE (following approach of
Philip Rosenau and Eli Turkel, Physica Scripta Vol. 31, 207-209,
1985). Unfortunatelly so far I had no expierence dealing with PDE
solvers. My attempts follows
(* initialization part *)
vp=ViewPoint/.Options[Plot3D];
fitparametersBaPP={a->0.9636882463878201,b->-99.380833179029,c->4.728158060984305};
fitparametersBaP={a->0.9920813114501638,b->-72.49635350172733,c->3.3297289189080654};fitCurvesfor3D[t_?NumberQ]:=Graphics3D[{fitCurves[[1,1]]/.{x_?NumberQ,y_?NumberQ}:>{t,x,y}}]
Show[{fitCurves=Plot[{a*Exp[b*x^2+c*x]/.fitparametersBaPP,a*Exp[b*x^2+c*x]/.fitparametersBaP},{x,-0.05,0.2},PlotStyle->{{Red,Thick},{Blue,Thick}}]}]
(* nonlinear PDE *)
equationsForSolution={
(-D[(diffusionC[1]*(((1 + \[Alpha][1])*n[1][t, x] + 2*n[2][t,
x])*Derivative[0, 1][n[1]][
t, x] + 2*\[Alpha][1]*n[1][t, x]*Derivative[0, 1][n[2]][t, x]))/
(n[1][t, x] + 2*n[2][t, x]), x] + Derivative[1, 0][n[1]][t,
x])/n0[1] == 0,
(-D[(diffusionC[2]*(2*\[Alpha][2]*n[2][t, x]*Derivative[0, 1][n[1]][t, x] +
(n[1][t, x] + 2*(1 + 2*\[Alpha][2])*n[2][t,
x])*Derivative[0, 1][n[2]][t, x]))/
(n[1][t, x] + 2*n[2][t, x]), x] + Derivative[1, 0][n[2]][t,
x])/n0[2] == 0};
Notations: n[ ][t, x] are ion (2 species) concentrations we are
interesting in, \[Alpha][ ] are some temperature ratio (constant close
to 1) , diffusionC[ ] are diffusion coefficients ( not well known, we
assume it should be between 100-500). The species initial
concentration ratio is about 100. We assume the following parameter
values
diffusionC[1]=100;diffusionC[2]=100;
n0[1]=100;n0[2]=1;
\[Alpha][1]=1;\[Alpha][2]=1;
bp=6/10 (*cm*);
tend=1/3000(*s*);
We start to solve from experimentally measured concentration
distribution (say, at t=0) and assume that at x=bp = 6/10 boundary
point all particles are absorbed.
initialCondition["GaussFit",1]=( n[1][0, x] ==
(n0[1]/diffusionC[1])*(a*Exp[b*x^2+c*x]/.fitparametersBaP));initialCondition["GaussFit",2]=(n[2][0, x] ==
(n0[2]/diffusionC[2])*(a*Exp[b*x^2+c*x]/.fitparametersBaPP));
As a result we have a solution
sol["GaussFit1"]=NDSolve[Flatten[{equationsForSolution,
initialCondition["GaussFit",1],initialCondition["GaussFit",2],n[1][t,bp] == 0,
n[1][t,-bp] == 0,n[2][t,bp] == 0,n[2][t,-bp] ==
0}],{n[1],n[2]},{t,0,tend},{x,-bp,bp},Method->{"MethodOfLines",
"SpatialDiscretization"->{"TensorProductGrid","DifferenceOrder"->"Pseudospectral"}}]
sprend["GaussFit1"]=GraphicsGrid[{{Plot3D[First[n[2][t,x] /.
sol["GaussFit1"]],{t,0,tend},{x,-bp, bp},
PlotRange->All,AxesLabel->{"t","x"},ViewPoint->Dynamic[vp]],Plot3D[First[n[1][t,x]/. sol["GaussFit1"]],{t,0,tend},{x,-bp, bp},
PlotRange->All,AxesLabel->{"t","x"},ViewPoint->Dynamic[vp]]}},ImageSize->1000]
The problems start when we try to solve this equations backward in
time. We need this because fitted initial distribution was measured
with some delay. So we want to look how this distribution might look
fraction of second before. Unfortunatelly Mathematica PDE solver refuses to
solve that way, claiming that singularity had formed. Sure it have to
form at some moment, however time scale looks unrealistic (Mathematica solves
backward to 10^-6 s only). So, my questions are
1) Does it is possible to get some solution, say 1/1000 fraction of
second before?
2) We have small boundary and initial condition mismatch at t=0 (the
function is not strictly 0 at point x=bp). How similar mismaches
practically are resolved in other similar PDE problems?
3) And the naive one. Can this nonlinear problem be classified as
being "elliptic", "parabolic" or some other kind, or this
classification is not possible for general nonlinear PDEs. In other
words, can linearization in principle help to determine the kind of
this PDE? Does this would help to choose the right options in Mathematica?
Sincerely, Arturas Acus