# Ndsolve of random (periodic) initial conditions via BSplineFunction

Suppose we want to solve the heat equation for a random initial temperature field T(t,x,y). Therefore I define the initial conditions with the BSplineFunction which nicely incorporates periodic boundary conditions and try to solve it with NDSolve:

Clear["Global*"]
Lx = 1; Ly = 1;
eq = D[T[t, x, y], t] == Laplacian[T[t, x, y], {x, y}];
pbc1 = T[t, x, 0] == T[t, x, Ly];
pbc2 = T[t, 0, y] == T[t, Lx, y];
T0[x_, y_] =
BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True][x,
y]
Plot3D[T0[x, y], {x, 0, Lx}, {y, 0, Ly}, PlotRange -> All,
AxesLabel -> {"x", "y", "T"}, ColorFunction -> "DarkRainbow",
Mesh -> All, MeshStyle -> Opacity[.2]]
T0[0, 0]
ic = T[0, x, y] == T0[x, y]
ic2 = T[0, x, y] == 0.5*Exp[-4 ((x - Lx/2)^2 + (y - Ly/2)^2)];
Monitor[AbsoluteTiming[
sol = NDSolve[{eq, pbc1, pbc2, ic},
T, {t, 0, 1}, {x, 0, Lx}, {y, 0, Ly},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 25}},
EvaluationMonitor :> (currentTime =
Row[{"t = ", CForm[t]}])]], currentTime]


The code unfortunately throws an error "Data … is not a rectangular tensor with dimensions …". I assume this results from the fact that T0 evaluated e.g. at grid point [0,0] returns a list {} and NDSolve probably needs just the float value.

I tried T0[0,0][] which yields the expected value. But this does not work for the definition of the initial condition: ic = T[0, x, y] == T0[x, y][], where x,y are not set yet. This returns just x, in reality it should first put in the values of x,y and only after that take list element 1 (if I am not mistaken).

Apart from that the code works fine which can be verified if you simply use ic2 instead of ic as initial condition in NDSolve.

Can someone help to get this running?

• Does ic := T[0, x, y] == T0[x, y][] help? (with SetDelayed) May 17 at 23:27
• First, transform the BSpline into an InterpolationFunction by: FunctionInterpolation[ BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True][x, y], {x, 0, 1}, {y, 0, 1}][x, y]. Further, your initial condition T0[x,y] is incredible complicated. THis will give numerical problems. May 18 at 8:39

Thank you Daniel, your comment solved the problem:

Clear["Global*"]
Lx = 1; Ly = 1;
eq = D[T[t, x, y], t] == Laplacian[T[t, x, y], {x, y}];
pbc1 = T[t, x, 0] == T[t, x, Ly];
pbc2 = T[t, 0, y] == T[t, Lx, y];
T0[x_, y_] =
FunctionInterpolation[
BSplineFunction[RandomReal[1, {30, 30, 1}], SplineClosed -> True][
x, y], {x, 0, 1}, {y, 0, 1}][x, y];
Plot3D[T0[x, y], {x, 0, Lx}, {y, 0, Ly}, PlotRange -> All,
AxesLabel -> {"x", "y", "T"}, ColorFunction -> "DarkRainbow",
Mesh -> All, MeshStyle -> Opacity[.2]]
ic = T[0, x, y] == T0[x, y];
Monitor[AbsoluteTiming[
sol = NDSolve[{eq, pbc1, pbc2, ic},
T, {t, 0, 1}, {x, 0, Lx}, {y, 0, Ly},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MaxPoints" -> 25}},
EvaluationMonitor :> (currentTime =
Row[{"t = ", CForm[t]}])]], currentTime]
Grid@Partition[
Table[Plot3D[T[t, x, y] /. sol, {x, 0, Lx}, {y, 0, Ly},
AxesLabel -> {"x", "y", "T"}, PlotPoints -> 31,
PlotRange -> {Automatic, Automatic, {0, 1}},
PlotLegends -> Automatic, ColorFunction -> "DarkRainbow",
Mesh -> All, MeshStyle -> Opacity[.2],
PlotLabel ->
Row[{Style["t=", Italic], Style[t, Italic]}]], {t, {0, 0.0005,
0.002, 1}}], 4] 