Help needed for understanding NDSolve architecture -Vacuum Heating Simulation

Question

I am a beginner in mathematica and I am facing a problem that is way bigger than my knowledge..

I am tryng to simulate a semplificated version of Brunel Effect, wich is a Laser-Plasma interaction effect, with an electrostatic Dawson's model.

The theory is simple: the plasma, with a sharp-step density shape, is hit by a sinusoidal laser impulse, some electrons will be pulled outside of the plasma and reach the vacuum region in which they are accelerated and then they re-enter into the plasma with more energy than before.

Mathematically the problem is also quite simple: the differential equations wich governs the electrons motion are different inside and outside the plasma, but for the consistency of the fluid model the electrons trajectories should never overlap and for avoiding this i have to "switch" the label of the trajectories when overlap occur. In a few words, i have to insert some checks while integrating the equations that recognize when an overlap occur and then restart the integration with the velocities swapped from that point, i.e. i have to change the initial conditions of the 2 electrons and then continue the integration till the next overlap.

The problem is that check should act inside the NDSolve procedure and while it is doing the integration.

Here the idea of the code i had been able to develop till now: Initialitation of the Field:

ClearAll;
 W = 10^6;
 tau = 10*Pi/W;
 c = 299792458;
 Am = 0.8*W*c;
 Ed[t_]:=If[t <((2Pi)/W),Sin[Pi*t/(2*tau)]*Sin[Pi*t/(2*tau)]*10*Am*Sin[W*t], Am*Sin[W*t]];
Plot[Ed[t], {t, 0, (2 Pi/W)*2}]

Define some constants and functions:

t =.
it = 0;
imax = 100;
check = 0;
check2 = 0;
aux = 0;
x0max = 4;
tmax = (2*Pi/W)*7;
V[x0_,t_]:=D[Z[x0,t],t];
x[x0_,t_]:=x0+Z[x0,t];

Initialitation of the eq.:

ndsdata = First@NDSolve`ProcessEquations[{D[V[x0, t], t] == 
  If[x0 + Z[x0, t] > 0, -10*W^2*Z[x0, t] - Ed[t], 
  10*W^2*x0 - Ed[t]], V[x0, 0] == 0, Z[x0, 0] == 0},  Z[x0, t], {x0, 0, x0max}, 
  {t, 0, tmax},  Method -> "ExplicitRungeKutta", MaxSteps -> Infinity]

Where Ed[t] is a Sinusoidal field, W is the pulse frequency, x0 the initial position of the electron, Z[x0,t] the electron displacement (the position of the electron at time t is x=x0+Z[x0,t] and V[x0,t] is the time derivate of Z[x0,t], which is the velocity of the electron

Iteration of the Solver:

For[i = 0, i < imax, i++, {it = it + 1/100000;
   NDSolve`Iterate[ndsdata, it];
 auxsol = NDSolve`ProcessSolutions[ndsdata];
   For[j = 0, j < x0max, j++,
If[Evaluate[x[j, it] /. auxsol] > 
  Evaluate[x[j + 1, it] /. auxsol],
 ndsdata = NDSolve`Reinitialize[ndsdata, {V[j, it] == V[j + 1, it], Z[j, it] ==
      Z[j + 1, it]}]];If[it < ((2 Pi/W))*7, imax = imax + 1, i = imax]]

And the process the solution and Plot:

sol = NDSolve`ProcessSolutions[ndsdata];

this code does not work, but if i try to integrate without the swap i am able to plot the (wrong) trajectories...

the code gives me no error, but it seems that the following If cycle neither became true nor false:

If[Evaluate[x[j, it] /. auxsol]>Evaluate[x[j + 1, it] /. auxsol],
{check = 1, V[j, it] == aux,V[j, it] = V[j + 1, it], V[j + 1, it] = aux,ndsdata =NDSolve`Reinitialize[
ndsdata, {V[j, it] == V[j + 1, it], Z[j, it] == Z[j + 1, it]}, 
check2 = 1]}]

So basically the swap (or what i have tried to do as a swap) never happens.

What i obtain are simply the solution of the integration, that are physically uncorrect:

Plot[Table[Evaluate[x[x0, t] /. sol], {x0, 0, x0max}],{t,0,(2*Pi/W)*7},
PlotRange -> {{0, ((2*Pi/W))*7}, {-700, 300}}]

What i am asking is: How can work inside Ndsolve procedure and change initial conditions of the equations at a certain time, when an event (in this case the overlap) occurs?

Thank you in advance and sorry for the expositive mess..

Answer 1

Solution with ParametricNDSolve

A fairly straightforward solution to this question can be obtained by observing that the underlying differential equation involves no derivatives with respect to x0, which therefore can be treated as a parameter. (Note that "StiffnessSwitching" should be used instead of "ExplicitRungeKutta".)

sol2 = ParametricNDSolve[{D[V[t], t] == 
    If[x0 + Z[t] > 0, -10*W^2*Z[t] - Ed[t], 10*W^2*x0 - Ed[t]], V[t] == D[Z[t], t], 
    V[0] == 0, Z[0] == 0}, {Z, V}, {t, 0, tmax}, {x0}, Method -> "StiffnessSwitching"];
Plot[Evaluate[Table[x0 + Z[x0][t] /. sol2, {x0, 0, x0max}]], 
    {t, 0, tmax}, AxesLabel -> {"t", "x"}]

The question, however, requires that the lines not cross and proposes accomplishing this by swapping the values of the lines when they do. Here, the lines are sorted by

f = Table[x0 + Z[x0][t] /. sol2, {x0, 0, x0max}, {t, 0, tmax, tmax/100}];
g = Map[Sort, f // Transpose, {1}] // Transpose;
h[x0_] := Interpolation[{Table[t, {t, 0, tmax, tmax/100}], g[[x0 + 1, All]]} 
    //Transpose]
Plot[Evaluate[Table[h[x0][t], {x0, 0, x0max}]], {t, 0, tmax}, 
    AxesLabel -> {"t", "x"}]

Higher resolution in x0 can be achieved in a straightforward manner, if desired.

Solution with NDSolve Components

The question proposed using NDSolve Components to obtain a solution. Doing so is not at all straightforward, in my view, but is interesting. One approach is

res = {};
ndsdata = First@NDSolve`ProcessEquations[{D[V[x0, t], t] == 
    If[x0 + Z[x0, t] > 0, -10*W^2*Z[x0, t] - Ed[t], 10*W^2*x0 - Ed[t]], 
    V[x0, t] == D[Z[x0, t], t], V[x0, 0] == 0, Z[x0, 0] == 0}, {Z, V}, {x0, 0, x0max}, 
    {t, 0, tmax}, Method -> "StiffnessSwitching"];
Do[it = i tmax/imax; NDSolve`Iterate[ndsdata, it]; 
    sol = NDSolve`ProcessSolutions[ndsdata]; 
    resi = Table[x[x0, it] /. sol, {x0, 0, x0max}]; 
    res = Append[res, resi]; ord = Ordering[resi];
    If[ord != Range[1, x0max + 1], 
      sav = Table[{Z[j, it], V[j, it]} /. sol, {j, 0, x0max}][[ord]];
      zf = Interpolation[sav[[All, 1]]]; vf = Interpolation[sav[[All, 2]]];
      ndsdata = First@NDSolve`Reinitialize[ndsdata, 
      {Z[x0, it] == zf[x0 + 1], V[x0, it] == vf[x0 + 1]}]], {i, imax}]
ListPlot[res // Transpose, DataRange -> {0, tmax}, Joined -> True, 
    AxesLabel -> {"t", "x"}]

Several points are worth making:

• This plot differs modestly from the previous plot at larger t, due to a loss of accuracy caused by 14 calls to Reinitialize with resorted initial conditions. The calculation takes about six seconds. Calling Reinitialize at every time step adds neither to the inaccuracy nor to the run-time.

• V[x0_,t_]:=D[Z[x0,t],t];, defined in the question, must be discarded, and an equivalent equation included in the argument of ProcessEquations.

• Dependent variables should be specified as {Z, V} rather than Z[x0, t]. The latter choice prevented the If in the question from executing.

• Even though the documentation might seem to suggest that

  ndsdata = NDSolve`Reinitialize[...]
  

  is correct, in fact

  ndsdata = First@NDSolve`Reinitialize[...]
  

  is necessary.

• Reinitialize followed by Iterate does not append additional solution results to existing results but instead discards the former. Hence, the earlier solution must be saved before calling Reinitialize.

• Specifying new initial conditions only at discrete values of x0 before calling Reinitialize, causes the dependent variables to be reinitialized to zero. Interpolation functions are used to here to provide continuous values to the reinitialized dependent variables.
• Consistent with the choice made in the question, x0 is sampled only at integer values during reinitialization. In practice, much higher resolution would be needed. Presumably, run-time increases linearly.

Solution with NDSolve WhenEvent

As suggested by user21, WhenEvent also might be used to trigger sorting of streamlines, although I have not attempted to do so.