For some reasons I need to implement the method of lines instead of using the NDSolve function directly. As described in this tutorial, this is how the simplest PDE can be solved:

H = 1.0; n = 10; dx = H/n;

U[t_] = Table[Subscript[u, i][t], {i, 0, n}];

eqns = Thread[
   D[U[t], t] == 
    Join[{D[Sin[2 π t], t] + (Sin[2 π t] - Subscript[u, 0][t])}, 
     ListCorrelate[{1, -2, 1}/dx^2, U[t], {1, 2}, {Subscript[u, n - 1][t]}]/8]];

initc = Thread[U[0] == Table[0, {n + 1}]];

lines = NDSolveValue[{eqns, initc}, U[t], {t, 0, 4}];

ParametricPlot3D[Evaluate[Table[{i dx, t, lines[[i + 1]]}, {i, 0, n}]], {t, 0, 4}, 
 PlotRange -> All, AxesLabel -> {"x", "t", "u"}]

A a result, functions $u_i(t)$ (i.e. "lines") in certain "space" points are obtained. NDSolve gives these functions as InterpolatingFunction objects. The solution in between lines can be found by interpolation.

How to get a two-dimensional InterpolatingFunction object from the functions $u_i(t)$, similarly to what NDSolve gives if I solve the equation directly?

I can create an array of point from the InterpolatingFunction and then interpolate, as follows

array = 
  Table[{{i dx, t1}, lines[[i + 1]] /. t -> t1}, {i, 0, n}, {t1, 0, 4, 0.1}];

function = ListInterpolation[array];

Plot3D[function[x, t], {t, 0, 4}, {x, 0, H}, PlotRange -> All]

but is there a better options? Can I somehow use the InterpolatingFunction object directly without introducing this additional division into points? Is there a build-in function to do this? I think my solution is not very effective because at first NDSolve interpolates between calculation points (when giving ODEs solutions), that I take this interpolation, and divide it into points again, and then, again, I interpolate...

  • 2
    $\begingroup$ Have a look at the tutorial section Transient PDEs. There you'll find a slightly different (but efficient) approach that may suite your needs. $\endgroup$
    – user21
    Jun 8, 2017 at 16:53

1 Answer 1


If the FEM method suggested by @user21 does not work, here is one way to interpolate between the line solutions. Since you've got lines in terms of the global variable t, we have to take some care with it (e.g., the use of Block[{t = t0}... in sol2d). It does an order-3 interpolation between the values of the four lines solutions whose x-coordinates are nearest the argument x.

xx = Range[0, n] dx;
xnf = Nearest[xx -> Automatic];

sol2d[x_?NumericQ, t0_?NumericQ] := Block[{t = t0},
  InterpolatingPolynomial[Transpose[{xx[[#]], lines[[#]]} &@xnf[x, 4]], x]]

Plot3D[sol2d[x, t], {x, 0, 1}, {t, 0, 4}, PlotRange -> All, 
 Mesh -> {xx}, MeshFunctions -> {#1 &}, MaxRecursion -> 3, 
 BoxRatios -> Automatic]

Mathematica graphics

It's about 3 times slower than the OP's ListInterpolation, but it uses the values from lines at time t instead of re-interpolating. This is slightly faster (~5%), since it avoids evaluating all the lines at t0:

sol2d[x_?NumericQ, t0_?NumericQ] := Block[{t},
  With[{xi = xnf[x, 4]},
   With[{yy = lines[[xi]]},
    t = t0;
    InterpolatingPolynomial[Transpose[{xx[[xi]], yy}], x]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.