# Help needed for understanding NDSolve architecture -Vacuum Heating Simulation

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@NDSolveProcessEquations[{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;
NDSolveIterate[ndsdata, it];
auxsol = NDSolveProcessSolutions[ndsdata];
For[j = 0, j < x0max, j++,
If[Evaluate[x[j, it] /. auxsol] >
Evaluate[x[j + 1, it] /. auxsol],
ndsdata = NDSolveReinitialize[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 = NDSolveProcessSolutions[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 =NDSolveReinitialize[
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..

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Sep 7 '15 at 17:04
• Have a look at the WhenEvent documentation. – user21 Sep 7 '15 at 18:14
• Your code is missing several constants, which the reader will need to attempt to reproduce and solve your problem. Also, what is the equation for Z[x0, t]? Finally, it what sense doe your code not work? If you are receiving error messages, what is the first such message? – bbgodfrey Sep 7 '15 at 19:46
• Hi, and thank you for the answer, i'll try to explain: I have reported only the code part concerned with the NDSolve Procedure, before I only Initialize the field Ed[t], some constans and then I define x[x0,t]=x0+Z[x0,t] as the spatial variation of the electron with initial position x0 and V[x0,t] as the partial derivative respect to time of Z[x0,t], which is the velocity. Z[x0,t] is the function I want to obtain with NDSolve, the differential equation which descrive its evolution is the one inside NDSolve: 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] – Cosimo Sep 7 '15 at 20:51
• Without the constants, expression for Ed[t], and equation for Z, it is unlikely that you will receive useful help. x[x0,t]=x0+Z[x0,t] does not define Z, because x also is undefined. Please add all this to your question. By the way, use @bbgodfrey in any comments you wish me to see promptly. I should add that you may be making this problem too hard for yourself. – bbgodfrey Sep 7 '15 at 21:15

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, Z == 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@NDSolveProcessEquations[{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; NDSolveIterate[ndsdata, it];
sol = NDSolveProcessSolutions[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@NDSolveReinitialize[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 = NDSolveReinitialize[...]


is correct, in fact

ndsdata = First@NDSolveReinitialize[...]


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.

• ok, I am impressed. Your answer is very deep and articulated, and will be a subject of study for me! Thank you for all your efforts, you really helped me in understaing this software! – Cosimo Sep 12 '15 at 15:27
• i have a doubt: i can see that there is no change in the trajectories between the swapped and unswapped ones..but there should be some difference, in fact the swaps concerns only with the velocities of the particles when they met each other ( i mean, when x0+Z[x0,t]=x1+Z[x1,t] ,for example, i have V[x0,t]=V[x1,t] and vice versa, but since the diff. eq. depends on x0, the evolution of the particle should be different from the one started from x1. In few words, we need to swap velocities, not the entire path.. Sorry if i hadn't been clear before – Cosimo Sep 12 '15 at 16:46
• @Cosimo I am no longer certain that I understand what you are seeking. It is clear that you wish to swap V, but to you also wish to swap Z, to do nothing else, or still some third action? – bbgodfrey Sep 12 '15 at 22:40
• now that i see the results, probably i was wrong in swapping Z, instead the only important swap may be V, keeping in memory the starting x0 even after the swap, thats why with the swap the results should be physically correct.. I am very sorry for the mess, i am still studing the thoery behind this problem while trying to perform this simulations! btw you really helped me so much, thank you again! – Cosimo Sep 13 '15 at 18:53
• @Cosimo I shall think about this some too, but not for a few days. I am traveling at the moment. Please let me know what you determine. Best wishes. – bbgodfrey Sep 14 '15 at 2:27