How to interpolate between InterpolatingFunction objects

Question

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...

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]

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]
    ]]]