NDSolve: Normalizing at every step

Question

Suppose I have an transport equation with an initial conditions:

sol = NDSolve[{D[y[x, t], t] - 4 D[y[x, t], x] == 0, 
              y[x, 0] == 1/Sqrt[2 \[Pi]] Exp[-(x)^2/2], 
              y[10, t] == 0}, 
              y[x, t], {x, -10, 10}, {t, 0, 5}, 
              MaxStepSize -> 0.1];
Plot3D[Evaluate[y[x, t] /. sol], {x, -10, 10}, {t, 0, 2.5}, 
       PlotRange -> All]

Suppose the packet is actually a probability distribution. The packet is moving to the left, and at some later time, part of the probability density goes out of the left boundary such that the sum of the probability density in the domain [-10,10] is less than 1.

My question is, what should I do such that, at every time step, the sum of the probability density in the domain [-10,10] is 1?

Answer 1

Here is the solution :

CRK4modified[]["Step"[rhs_, t_, h_, y_, yp_]] := 
 Module[{k0, k1, k2, k3, ampli},
  k0 = h yp;
  k1 = h rhs[t + h/2, y + k0/2];
  k2 = h rhs[t + h/2, y + k1/2];
  k3 = h rhs[t + h, y + k2];
  ampli = 
   If[t < 0.2, 1,  6.5 (* 6.5 : manually trimmed ! *)/Total[y]];
  {h, -y + ampli ((k0 + 2  k1 + 2  k2 + k3)/6 + y)}]
CRK4modified[___]["DifferenceOrder"] := 4
CRK4modified[___]["StepMode"] := Fixed

tEnd = 3;
sol20 = NDSolve[
         {  D[y[x, t], t] - 4 D[y[x, t], x] == 0, 
            y[x, 0] == 1/Sqrt[2 \[Pi]] Exp[-x^2/2],
            y[10, t] == 0}, 
         y,
         {x, -10, 10},
         {t, 0, tEnd},
         Method -> {MethodOfLines, 
            Method -> {"FixedStep", "StepSize" -> 0.02,Method -> CRK4modified},
            TemporalVariable -> t
            }
         ][[1]]

Plot3D[Evaluate[y[x, t] /. sol20], {x, -10, 10}, {t, 0, tEnd}, 
 PlotRange -> {Automatic, Automatic, {-0.1, 1}}]

The solution is based on the use of a do-it-yourself algorithm for solving the ODEs used by the PDE. I have token the example of the Runge-Kutta algorithm (CRK4) from the documentation (here, mathematica version 8) and modified it by dividing the solution-vector y (which values are > 0) by Total[y]. This maintains constant the area under the curve at each step.

This algorithm is intentionally effective only after t=0.2 second. This permits to adjust a multiplicative constant coefficient on Total[y] (6.5 here, try 8 to see the effectiveness of the algorithm, you will see that the height of the wave has a step at t=0.2)

Few details about the code of the modification of the Runge-Kutta algorithm CRK4 :

We want to divide the solution-vector by the surface of the curve. The correction factor is ampli = If[t < 0.2, 1, 6.5/Total[y]]; Total[y] is proportional to the surface, with a undetermined constant value that I have manually trimmed (the 6.5, which depends on the number of points of the spatial grid)

The term -y + ampli ((k0 + 2 k1 + 2 k2 + k3)/6 + y) in the increment of the solution-vector. As we want to multiply the vector (not the increment) by ampli, starting from the normal CRK4 increment (k0 + 2 k1 + 2 k2 + k3)/6, I have done a translation +y -> multiply by ampli-> back-translation (-y).

Answer 2

NOT a answer.
Here is a no successful attempt. Definitive ?

Mathematica has utilities that permit the user to manage time during temporal simulations. One can advance forward (and even backward) on a certain amount of time. Between theses intervals, I hope it should be possible to change the data.

It is documented in http://reference.wolfram.com/mathematica/tutorial/NDSolveStateData.html#93858351.

I currently use Mathematica 8.0.4., but the functionality I use here are old (maybe Version 5). It is possible that I miss a more recent functionality.

Here is a prolog to this simulation :

ndssdata = First[NDSolve`ProcessEquations[{
     D[y[x, t], t] - 4 D[y[x, t], x] == 0,
     y[x, 0] == 1/Sqrt[2 \[Pi]] Exp[-(x)^2/2],
     y[10, t] == 0
     },
   y,
   {x, -10, 10},
   {t, 0, 100}, MaxStepSize -> 0.01, 
   Method -> {MethodOfLines, TemporalVariable -> t}]]

then one can simulate the first 1.1 seconds :

NDSolve`Iterate[ndssdata, 1.1]

here is the code to see the result :

ndsol = NDSolve`ProcessSolutions[ndssdata]
Plot3D[Evaluate[y[x, t] /. ndsol], {x, -10, 10}, {t, 0, 1.1}, 
 PlotRange -> All]

We can get the actual data from inside the simulator :

actualData=ndssdata@"SolutionVector"["Forward"];
ListPlot[actualData]

It seems that we can modify them, for example multiply them by 2 :

    Unprotect[NDSolve`StateData]
    ndssdata@"SolutionVector"["Forward"]= 2 ndssdata@"SolutionVector"["Forward"]
    ndssdata@"SolutionDerivativeVector"["Forward"]= 2 ndssdata@"SolutionDerivativeVector"
["Forward"]
    Protect[NDSolve`StateData]

but it is not effective : ndssdata@"SolutionVector"["Forward"] is not modified

EDIT

At the moment, it seems that to write in ndssdata@"SolutionVector"["Forward"]is not a solution.

If one want to stop the simulation, change the data, and resume the simulation, there is the possibility to use NDSolve``ReInitialize :

newstate = 
 First @ NDSolve`Reinitialize[
   ndssdata, {y[x, 0] == 2((ndsol[[1, 2]])[x, 2])}] (* 2="normalisation"*)

NDSolve`Iterate[newstate, 2]

Result :

ndsolnew = NDSolve`ProcessSolutions[newstate]
Plot3D[Evaluate[y[x, t] /. ndsolnew], {x, -10, 10}, {t, 0, 2}, 
 PlotRange -> All]

we see the continuation of the simulation.

Not very good because :
- the time reset to 0 at every ReInitilization
- not applicable to each step

We can view the junction of the two parts with the following code (merge the whole + time of first simulation elongated to 2 seconds + no factor 2):

(* part 1 : *)
ndssdata = 
  First[NDSolve`ProcessEquations[  
      {  D[y[x, t], t] - 4 D[y[x, t], x] == 0,
         y[x, 0] == 1/Sqrt[2 \[Pi]] Exp[-(x)^2/2],
         y[10, t] == 0},  
      y,
      {x, -11, 10},
      {t, 0, 100},
      MaxStepSize -> 0.01,
      Method -> {MethodOfLines, TemporalVariable -> t}
      ]];
NDSolve`Iterate[ndssdata, 2]
ndsol = NDSolve`ProcessSolutions[ndssdata]
gr1 = Plot3D[Evaluate[y[x, t] /. ndsol],
                {x, -10, 3},{t, 0, 2}, 
                PlotRange -> All,PlotStyle -> Specularity[White, 5]];
tEndPartOne = Last[ndssdata @ "CurrentTime"];  

(* part 2 : *)
newstate = 
 First@NDSolve`Reinitialize[
   ndssdata,  
   {y[x, tEndPartOne] == ((ndsol[[1, 2]])[x, 2])}
   ] 
NDSolve`Iterate[newstate, 3.7]
ndsolnew = NDSolve`ProcessSolutions[newstate]
gr2 = Plot3D[Evaluate[y[x, t] /. ndsolnew],
             {x, -10, 3},{t, 2, 3.7}, 
             PlotRange -> All];  

(* whole :*)
Show[gr1, gr2]