# How to use NDSolve with limited RAM?

How can I use NDSolve when running out of RAM?

The code I am implementing:

sol = NDSolve[ pdes, {A, B, C},
{x, -L, L}, {t, Tf, Tf} ]


in order to obtain the end solution of NDSolve to be exported to the next run (due to limited RAM).

In previous attempts the solution of NDSolve was exported as an interpolating function, such as:

Export["sol_Tf.wdx",sol]


sol_Tf.wdx is the export file of the NDSolve run.

Also see my previous question.

I then imported it as:

Import["sol_Tf.wdx"]


and used the boundary conditions:

A[x, Ti] == A[x, Tf]
B[x, Ti] == B[x, Tf]
C[x, Ti] == C[x, Tf]


But the results were different from running it all in one go. What went wrong here? Is there a better way of doing this?

For example, in the heat equation:

sol = NDSolve[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0,
u[t, 0] == Sin[t], u[t, 5] == 0}, u, {t, 0, 10}, {x, 0, 5}]


I would want to continue it with the end results of the previous run, as the new boundary conditions of the next run; RAM free, and by running it again from t = 10 to t = 20.

• I think it would be useful if you added a simple PDE exmaple such that people can experiment with this. As the question stands now, I'd first need to come up with an example to try this. – user21 Sep 16 at 5:56

You may use the NDSolve Components and Data Structures tutorial to control the memory usage and save down the state of intermediate runs.

Initialise the NDSolveStateData for the complete range you need to create solutions. Below is done for 0 <= t <= 30 and I will iterate in chunks of 10.

ndsStateData = First@NDSolveProcessEquations[
{
D[u[t, x], t] == D[u[t, x], x, x],
u[0, x] == 0,
u[t, 0] == Sin[t],
u[t, 5] == 0
},
u,
{t, 0, 30}, {x, 0, 5}
];


Next I iterate in 3 chunks of 10. ndsStateData can be saved down after the current chunk's solution is extracted with NDSolveProcessSolutions. Below I reassign ndsStateData with the reinitialised NDSolveStateData instead.

sols = {};
Module[{step = #},
NDSolveIterate[ndsStateData, {(step - 1)*10, step 10}];
AppendTo[sols, u /. NDSolveProcessSolutions[ndsStateData]];
ndsStateData =
First@NDSolveReinitialize[ndsStateData, {u[step 10, x] == sols[[step]][step 10, x]}];
] & /@ Range@3;


sols contains the 3 solutions over each chunk. Notice the difference in the domains.

sols Plotting the solutions.

Plot3D[sols[[#]][t, x], {t, (# - 1) 10, # 10}, {x, 0, 5},
PlotRange -> Full,
PlotStyle -> ColorData[#]
] & /@ Range@3 // Show Hope this helps.

The following throws some inconsistency warnings, but appears to work (i.e., no obvious discontinuities at the seam). Also, I did not change the boundary condition, but I fed the final condition of the first solution to the "initial" condition of the second solution.

uif1 = NDSolveValue[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0,
u[t, 0] == Sin[t], u[t, 5] == 0}, u, {t, 0, 10}, {x, 0, 5}];
uif2 = NDSolveValue[{D[u[t, x], t] == D[u[t, x], x, x],
u[10, x] == uif1[10, x], u[t, 0] == Sin[t], u[t, 5] == 0},
u, {t, 10, 20}, {x, 0, 5}];
plt1 = Plot3D[uif1[t, x], {t, 0., 10.}, {x, 0., 5.},
PlotRange -> {{0, 20}, {0, 5}, {-1, 1}}];
plt2 = Plot3D[uif2[t, x], {t, 10., 20.}, {x, 0., 5.},
PlotRange -> {{0, 20}, {0, 5}, {-1, 1}},
PlotStyle -> RGBColor[0.22, 1., 0.44]];
Show[{plt1, plt2}] 