1
$\begingroup$

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][[1]] >= 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.

I have tried to add the additional boundary conditions:

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – yohbs Feb 2 '15 at 13:36
  • $\begingroup$ You have mismatched quotation marks around MinPoints in your 2-D version. $\endgroup$ – bbgodfrey Feb 2 '15 at 13:48
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Feb 2 '15 at 14:53
  • $\begingroup$ @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. $\endgroup$ – Enter Feb 2 '15 at 15:18
  • $\begingroup$ @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! $\endgroup$ – Enter Feb 2 '15 at 15:55
1
$\begingroup$

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}][[1]] >= 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]

enter image description here

$\endgroup$
  • $\begingroup$ I just ran the case without difficulty. $\endgroup$ – bbgodfrey Feb 3 '15 at 4:01
  • $\begingroup$ @ 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}][[2]]; ymax = y /. Maximize[{h[x, y, t], 0<x<L && 0<y<L && 0<t<= tmax}, {x, y, t}][[2]];.Rewriting WhenEvent as WhenEvent[h[xmax,ymax,t]>=p, "StopIntegration"]. But I fail to combine them. $\endgroup$ – Enter Feb 3 '15 at 4:08
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Feb 3 '15 at 5:17
  • $\begingroup$ @ bbgodfrey, Thanks again. I will handle it and I am using Mathematica 9.0. $\endgroup$ – Enter Feb 3 '15 at 5:46
  • $\begingroup$ Works equally well for me with Mathematica 9.0.1.0, although the calculation stops slightly earlier, at 119.10119006911476. $\endgroup$ – bbgodfrey Feb 3 '15 at 6:08

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.