# NDSolveReinitialize start from previous solution but with different boundary conditions

I am solving a set of coupled partial differential equations (some variation of a 1 dimensional compressible fluid flow). I want to solve the equations (NDSolve, this works fine) but then continue solving with different boundary conditions. I am trying to use NDSolveReinitialize for this but it gives me an error.

Current code: I first have someFunction as boundary condition. Then I want to solve the system, and continue with NEWFUNCTION as boundary condition.

(* rhostart, vstart, estart are defined previously *)
boundaryconditions = (
rho[0, z] == rhostart[z]
&& v[0, z] == vstart[z]
&& e[0, z] == estart[z]
&& ...
&& e[t, 0] == someFunction[t]
&& ...);

statedata = First[NDSolveProcessEquations[
{eq1, eq2, eq3, boundaryconditions},
{rho, v, e}, {t, 0, tmax}, {z, 0, zmax}
]]

(* I call Iterate in a loop to check if the result is sufficiently stable (I omitted that code here) *)
NDSolveIterate[statedata, curTime]

(* When stable, I process the result *)
solution = NDSolveProcessSolutions[statedata]
rhostart[z_] := Evaluate[rho[curTime,z] /. solution];
estart[z_] := Evaluate[e[curTime,z] /. solution];
vstart[z_] := Evaluate[v[curTime,z] /. solution];

(* I compute a few things using the result, omitted here *)

(* Now I want to change a boundary condition by changing e[t,0] *)
newboundaryconditions = (
rho[0, z] == rhostart[z]
&& v[0, z] == vstart[z]
&& e[0, z] == estart[z]
&& ...
&& e[t, 0] == NEWFUNCTION[t]
&& ...);
NDSolveReinitialize[statedata, {newboundaryconditions}];

This gives me the following error:

NDSolveReinitialize::ndsv: Cannot find starting value for the variable v.

More details: One of my boundary conditions involves setting a variable at x=0 to a certain value. I then let the system evolve in time untill it reaches a steady state. Subsequently I want to compute a few values and then change this boundary condition to a different value and let the system evolve again in time starting from where it left off untill it reaches a sufficiently steady state again.

Important: I could change the original boundary condition function someFunction to incorporate all changes at long time intervals. However I do not want to keep the full solution of NDSolve in memory all the time (it already reaches 6 GB quickly) so I think it is neccesary to restart the NDSolve process each time.

• Are the new boundary conditions consistent with the previous solution? For example, if you had some change like f[0, x] == 0 to f[0, x] == 1, then NDSolve won't be able to reconcile the discontinuity. – 2012rcampion Feb 20 '15 at 14:10
• Yes. The initial function someFunction[t] interpolates from value A to B over time. The solution starts with e[0,0] = A and at the end of NDSolve, the solution has e[tfinal, 0] = B. The new boundary condition NEWFUNCTION is similar but interpolates from B to C, so it should be consistent. – Tom Bannink Feb 20 '15 at 15:16
• I was thinking the problem might be in the way that I set the new rhostart , vstart and estart objects since they are InterpolatingFunction objects so maybe I need to define them like rhostart[x_?NumericQ]:=... or something. I tried with and without using Evaluate but no luck so far. – Tom Bannink Feb 20 '15 at 15:27
• Try using Set instead of SetDelayed (e.g. vstart[z_] = v[curTime,z] /. solution). Also, maybe try saving a copy of statedata before iterating the solution, and using that in the Reinitialize. – 2012rcampion Feb 20 '15 at 15:36
• Thanks for the suggestions! I tried both suggestions (and different combinations of them) but sadly I still get the same error. I did find out the following: If I call NDSolve Reinitialize with the original boundary conditions as second argument, it gives the same error. I should have checked this before. The problem is probably not caused by the way I set the Interpolating function objects and so on. Why would it fail when I call Reinitialize with the same boundary conditions as in the original call? – Tom Bannink Feb 22 '15 at 21:41

I can confirm this problem appears when altering the boundary conditions, but not when altering the initial conditions. Here's an example of the former (the latter case is handled successfully in the documentation):

(* setup simple damped wave equation *)
diffeq = Derivative[0, 2][u][x, t] + Derivative[0, 1][u][x, t] -
Derivative[2, 0][u][x, t] == 0;
l = 5; tmax = 5; (* set bounds *)
(* set boundary and initial conditions *)
boundary1 = (u[0, t] == 0 && u[l, t] == Sin[t]^2 && u[x, 0] == 0 &&
Derivative[0, 1][u][x, 0] == 0);
(* create NDSolve data object *)
ndssdata =
First@NDSolveProcessEquations[{diffeq, boundary1},
u, {x, 0, l}, {t, 0, tmax}];
(* integrate and get solution*)
NDSolveIterate[ndssdata, tmax];
soln1 = NDSolveProcessSolutions[ndssdata];
(* setup new initial conditions *)
u0 = FunctionInterpolation[u[x, tmax] /. soln1, {x, 0, l, 0.01}][x];
du0 = FunctionInterpolation[
Derivative[0, 1][u][x, tmax] /. soln1, {x, 0, l, 0.01}][x];
(* set initial conditions and new boundary conditions *)
boundary2 = (u[0, t] == 0 && u[l, t] == Sin[t + tmax]^2 &&
u[x, 0] == u0 && Derivative[0, 1][u][x, 0] == du0);
(* error on reinitialize *)
ndssdata2 = First@NDSolveReinitialize[ndssdata, boundary2]
(* NDSolveReinitialize::ndsv: Cannot find starting value for the variable u. *)

As a work-around, you can just re-process the equations. If your individual integrations take a long time (I suspect they do if they generate ~6GB of data), and your equations are not hideously complex, then this should add minimal overhead.

ndssdata2 =
First@NDSolveProcessEquations[{diffeq, boundary2},
u, {t, 0, tmax}, {x, 0, l}];
NDSolveIterate[ndssdata2, tmax];
soln2 = NDSolve`ProcessSolutions[ndssdata2];
(* solutions match! *)
Plot3D[If[t <= tmax, u[x, t] /. soln1, u[x, t - tmax] /. soln2], {x,
0, l}, {t, 0, 2 tmax}]

If you're doing it this way, and you don't need the extra control over the solving process, then it's a bit simpler just to use regular NDSolve:

soln1 = First@NDSolve[{diffeq, boundary1}, u, {x, 0, l}, {t, 0, tmax}];
(* set u0, du0, and boundary2 as before *)
soln2 = First@NDSolve[{diffeq, boundary2}, u, {x, 0, l}, {t, 0, tmax}];
solnfull = {u ->
FunctionInterpolation[
If[t <= tmax, u[x, t] /. soln1, u[x, t - tmax] /. soln2], {x, 0,
l}, {t, 0, 2 tmax}]};
Plot3D[u[x, t] /. solnfull, {x, 0, l}, {t, 0, 2 tmax}]
• Thank you. It is indeed easier to simply re-process the equations. The processing does take a minute or two but this is still a small amount of time compared to the integrations which take an hour or so. I still use the Iterate method (as opposed to the regular NDSolve) because I want to be able to stop the process every now and then to check the intermediate results. – Tom Bannink Feb 25 '15 at 14:45