# NDSolve with initial condition: initial conditions did not evaluate to an array of numbers of depth 1(for 1D) or 2 (for 2D) on spatial grid

I have the following naive code to solve a PDE in two spatial dimensions (x,y) with periodic boundary conditions:

L = 20;
p = 8;
tmax = 500;
β = 1;
σ = 2;
myfun = First[h /. NDSolve[{D[h[x, y, t], t] +
Div[h[x, y, t]^3*Grad[Laplacian[h[x, y, t], {x, y}], {x, y}], {x, y}] +
Div[h[x, y, t]^3*Grad[h[x, y, t], {x, y}], {x, y}] -
Div[(h[x, y, t]^2*Grad[h[x, y, t], {x, y}])/(1 + β*h[x, y, t])^2, {x, y}] == 0,
h[x, y, 0] == 1 +
1/(2*π*σ^2)*Exp[-((x - 10)^2/(2*σ^2) + (y - 10)^2/(2*σ^2))],
h[0, y, t] == h[L, y, t], h[x, 0, t] == h[x, L, t],
WhenEvent[
NMaximize[{myfun[a, b, c], a >= 0 && a <= L, b >= 0 && b <= L, c > 0 && c <= t},
{a, b, c}, WorkingPrecision -> 10][] >= p, "StopIntegration"]},
h, {x, 0, L}, {y, 0, L}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 40, "MaxPoints" -> 40,
"DifferenceOrder" -> 4}},
AccuracyGoal -> 20, WorkingPrecision -> MachinePrecision, StepMonitor :> Print[t]]]


But when I try to run it, Mathematica prompts:

The initial conditions did not evaluate to an array of numbers of depth 2 on the spatial grid. Initial conditions for partial differential equations should be specified as scalar functions of the spatial variables.

I do not understand the first issue. However, as for the second issue I do think my initial condition

h[x, y, 0] == 1 + 1/(2*π*σ^2)*Exp[-((x - 10)^2/(2*σ^2) + (y - 10)^2/(2*σ^2))]


has be specified as scalar function of x, y.

Derivative[1, 0, 0][h][0, y, t] == Derivative[1, 0, 0][h][L, y, t],
Derivative[0, 1, 0][h][x, 0, t] == Derivative[0, 1, 0][h][x, L, t],


but the issues remain. Then I have tried to run a 1D version code:

myfun1D = First[h /.NDSolve[{D[h[x, t], t] + D[h[x, t]^3*D[h[x, t], {x, 3}], x] +
D[h[x, t]^3*D[h[x, t], x], x] -
D[(h[x, t]^2*D[h[x, t], x])/(1 + β*h[x, t])^2, x] == 0,
h[x, 0] == 1 + 1/(2*π*σ^2)*Exp[-((x - 10)^2/(2*σ^2))],
h[0, t] == h[L, t],
WhenEvent[
NMaximize[{myfun1D[a, c], a >= 0 && a <= L, c > 0 && c <= t}, {a, c},
WorkingPrecision -> 10][] >= p, "StopIntegration"]},
h, {x, 0, L}, {t, 0, tmax},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 40, "MaxPoints" -> 40,
"DifferenceOrder" -> 4}}, AccuracyGoal -> 20,
WorkingPrecision -> MachinePrecision, StepMonitor :> Print[t]]]


The error still persists.

Any ideas what may be causing these problems ? (If you see any other improvements I could make to my code I would be happy to hear them.)

Update: As suggested by @bbgodfrey, I changed WorkingPrecision->MachinePrecision. Now, it can work, but still with some warning about WhenEvent. To terminate the integration somehow, I employed WhenEvent to compare the max of myfun in space(x,y) and time(t) with a threshold p. Is there any other way to realize this stopping criterion?

• I would recommend to create a minimal example showing the problem. This way your pretty likely to find it yourself. For example, on my machine, if I delete the WhenEvent onwards, and work with the default options of NDSolve, the problem disappears. – yohbs Feb 2 '15 at 13:36
• You have mismatched quotation marks around MinPoints in your 2-D version. – bbgodfrey Feb 2 '15 at 13:48
• The error is associated with WorkingPrecision. Deleting or setting it to MachinePrecision eliminates the error. I find this very strange. Note that I deleted WhenEvent, because mdfun is undefined. – bbgodfrey Feb 2 '15 at 14:53
• @yohbs, Thanks a lot. I will try to simply my question. I posted the full problem because I was doubting I probably did something wrong elsewhere. Specially, I did not sure if I am using Div, Grad, and Laplacian properly in my 2D version. When I adopt your advice, it works. But I just want to explicitly control the method, grid, and accuracy or so. – Enter Feb 2 '15 at 15:18
• @bbgodfrey, thanks for your kind help! Yes, I missed a " before MinPoints. I adopted you advice just now, WorkingPrecision->MachinePrecision. It is very surprising. It does work, however, I use the similar accuracy level in my other 1D code. They all work properly. To terminate the integration somehow, I use WhenEvent to compare the max of myfun in space(x,y) and time(t) with a threshold p. It is a typo that I have written mdfun instead of myfun incorrectly.:) Is there any other way to realize this criterion? Thanks again! – Enter Feb 2 '15 at 15:55

As I mentioned in a Comment, the error is associated with WorkingPrecision. Deleting or setting it to MachinePrecision eliminates the error.

However, WhenEvent also is producing errors. I believe that it should be rewritten as

WhenEvent[NMaximize[{h[x, y, t]}, {x, y}][] >= p, "StopIntegration"]


However, even this does not work, perhaps because NDSolve is passing h to NMaximize without evaluating it as a function at t. So, I suggest that you instead find Max at several likely points, for instance,

WhenEvent[Max[h[10, 10, t], h[0, 0, t]] >= p, "StopIntegration"]


which works, terminating NDSolve at t = 119.1766609981875.

Plot3D[myfun[x, y, 119.17666], {x, 0, 20}, {y, 0, 20}, PlotRange -> All] • I just ran the case without difficulty. – bbgodfrey Feb 3 '15 at 4:01
• @ bbgodfrey,I tried it with WhenEvent[h[10, 10, t] >= p, "StopIntegration"]just a slight modification of yours with Max[h[10, 10, t], h[0, 0, t]] >= p. However, warning about WhenEvent remains:"The function value ...>=p is not True or False...". I am thinking is there any way to stop NDSolve like that:xmax=x/.Maximize[{h[x, y, t], 0<x<L && 0 < y < L && 0 < t <= tmax}, {x, y, t}][]; ymax = y /. Maximize[{h[x, y, t], 0<x<L && 0<y<L && 0<t<= tmax}, {x, y, t}][];.Rewriting WhenEvent as WhenEvent[h[xmax,ymax,t]>=p, "StopIntegration"]. But I fail to combine them. – Enter Feb 3 '15 at 4:08
• WhenEvent[h[10, 10, t] >= p, "StopIntegration"] works fine for me. Without seeing precisely what you are doing, I cannot offer further advice. However, I can say that that your alternative involving Maximize does not work. – bbgodfrey Feb 3 '15 at 5:17
• @ bbgodfrey, Thanks again. I will handle it and I am using Mathematica 9.0. – Enter Feb 3 '15 at 5:46
• Works equally well for me with Mathematica 9.0.1.0, although the calculation stops slightly earlier, at 119.10119006911476. – bbgodfrey Feb 3 '15 at 6:08