0
$\begingroup$

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][[1]] 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][[1]], 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?

$\endgroup$
2
  • $\begingroup$ Does ic := T[0, x, y] == T0[x, y][[1]] help? (with SetDelayed) $\endgroup$
    – MelaGo
    May 17 at 23:27
  • $\begingroup$ 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. $\endgroup$ May 18 at 8:39
1
$\begingroup$

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]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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