# Why is a bump function making NDSolve take forever to solve? [closed]

I am attempting to solve a system of relatively complicated elliptic PDEs via a relaxation technique. In particular, I am trying to solve $\textrm{div} LW = S$ where $S$ is some vector and $L$ is the operator $LW_{ij} = \nabla_i W_j + \nabla_j W_i -\frac{2}{n} \textrm{div}W g_{ij}$, essentially on the flat torus.

There's a bit more complication here, but there's a function I put in at the beginning that describes the geometry of the problem. Since it's on a flat torus, it needs to be a periodic function. Even fairly complicated trig functions like $\sin(x)\cos^3(y) + \sin (z)$ or something solve in 20 or 30 minutes.

However, when I put in a standard bump function, $\exp(\frac{-1}{1 - (x - 2)^2 - (y - 2)^2})$ in a circle around (2,2) and 0 elsewhere (so it's piecewise), it takes forever to solve. I left it open over the weekend and it still hasn't solved it.

Any suggestions?

For my code, the important bits (I think) are

vectorEquationStart =
Join[vectorEquation, {W1[time, 0, y] == W1[time, 2 π, y],
W1[time, x, 0] == W1[time, x, 2 π], W1[0, x, y] == Sin[x],
W2[0, x, y] == Sin[y], W3[0, x, y] == 0}];


which takes the equations I set up, vectorEquation, and attaches the boundary conditions I use. I am solving for W1, W2 and W3. Then

soln = NDSolve[
vectorEquationStart /. {α ->ξ}, {W1, W2, W3}, {time, 0, 10}, {x, 0, 2π}, {y, 0, 2π},
MaxStepFraction -> 1/25]


is what I try to solve. The MaxStepFraction is not necessarily the problem; I had another one running at the same time without it (slightly different problem, but similar) and it also has not given me an answer. The alpha in this code is a simple parameter (ξ=7/8 here). That has given be problems in previous problems if I make it too close to 1 because of how it affects the equations, but 7/8 has not given me problems before.

Is there something obvious I don't know or does anyone have any suggestions on how to make this run faster?

EDIT: I eventually got an answer from Mathematica, but it within 1.5 sec was up near $10^{50}$ or something. I snooped around and think that the problem may be how mathematica deals with my piecewise function near the boundary of the circle.

The function is mathematically smooth, but when I tried to find the min of $1-\frac{4 e^{\frac{2}{7-4x+x^2 - 4y +y^2}} ((x-2)^2+(y-2)^2)}{(7-4x+x^2-4y+y^2)^4}$, which should be positive everywhere for my calculations to work, I got $-1.5\cdot 10^{1123214}$ or so which could account for the observed behavior of the solver. If I look at the min on a slightly smaller circle, say $(x-2)^2+(y-2)^2 <.95$, it gives a reasonable answer, about .37.

Here's my code. Warning, it's about 40 lines to get down to where I actually run the NDSolve.

vecNorm[vec_?VectorQ, g_?MatrixQ] := Sqrt[{vec}.g.Transpose[{vec}]][[1, 1]]

Minkowski = {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0,0, -1}}

coord = {x, y, z}

f[x_, y_] := Piecewise[{{Exp[-1/(1 - (x - 2)^2 - (y - 2)^2)], (x - 2)^2 + (y - 2)^2 <
1}, {0, (x - 2)^2 + (y - 2)^2 >= 1}}];


t is the parametrization of the graph.

t[x_, y_, z_] := α f[x, y]


X, Y, Z are the three tangent directions

X = {1, 0, 0, D[t[x, y, z], x]}

Y = {0, 1, 0, D[t[x, y, z], y]}

Z = {0, 0, 1, D[t[x, y, z], z]}


A is the normal vector.

A = {D[t[x, y, z], x], D[t[x, y, z], y], D[t[x, y, z], z], 1}

Vect = {X, Y, Z, A}

NA = Simplify[vecNorm[A, -1*Minkowski]]


(The inside of the sqrt in NA MUST be positive or else I picked a bad function t(x,y,z) to start. It SHOULD work for this t(x,y,z), but Mathematica seems to mess up at the boundary of the piecewise parts.) H is the un-normalized 2nd fundamental form, calculated using Weingarter's rule.

H = -Table[Sum[D[A[[k]], coord[[k]]]*Vect[[i, k]]*Vect[[j, k]], {k, 3}], {i, 3}, {j, 3}]


h is the normalized 2nd fundamental form.

h = (1/NA)*H


g is the metric.

g = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}} - Table[D[t[x, y, z], coord[[i]]]*D[t[x, y, z], coord[[j]]], {i, 3}, {j, 3}]


G is the inverse metric.

G = Simplify[Inverse[g]]


τ is the trace of the second fundamental form, i.e. the mean curvature.

τ = Simplify[Tr[h.G]]


dτ is then the gradient of the mean curvature. We call that gradient T

T = Grad[τ, {x, y, z}];

affine := affine = Simplify[Table[(1/2)*Sum[(G[[i, s]])*(D[g[[s, j]], coord[[k]]] + D[g[[s, k]], coord[[j]]] - D[g[[j, k]], coord[[s]]]), {s, 1, 3}], {i, 1, 3}, {j, 1, 3}, {k, 1, 3}]]


This metric Divergence only works if the X, Y, Z are coordinates. (which they are, so no worries)

metricDiv[W_?VectorQ, g_?MatrixQ , coords_?VectorQ] := 1/Sqrt[Det[g]]*Sum[D[W[[i]] Sqrt[Det[g]], coords[[i]]], {i, 3}]

metricL[W_?VectorQ, g_?MatrixQ, Christoffel_, coords_?VectorQ] := Table[Sum[g[[m, j]]*(D[W[[m]], coords[[i]]] + Sum[Christoffel[[m, i, k]]*W[[k]], {k, 3}]), {m, 3}] + Sum[g[[i, m]]*(D[W[[m]], coords[[j]]] + Sum[Christoffel[[m, j, k]]*W[[k]], {k, 3}]), {m, 3}], {i, 3}, {j, 3}] - 2/3 metricDiv[W, g, coord]*g


The trace free part of h is given below as S.

S = Simplify[h - 1/3 g*τ];

DS = Table[metricDiv[S[[i]], g, coord], {i, 3}];


We now define W, LW, etc.

W = {W1[time, x, y], W2[time, x, y], W3[time, x, y]};

testW = {W1, W2, W3};

LW = metricL[W, g, affine, coord];

vectorEquation = Table[Derivative[1, 0, 0][testW[[i]]][time, x, y] ==  metricDiv[LW[[i]], g, coord] - DS[[i]], {i, 3}];

vectorEquationStart = Join[vectorEquation, {W1[time, 0, y] == W1[time, 2 π, y], W1[time, x, 0] == W1[time, x, 2 π], W1[0, x, y] == Sin[x], W2[0, x, y] == Sin[y], W3[0, x, y] == 0}];

We now attempt to solve the PDE we just set up. ξ will be α

ξ = 7/8;

soln = NDSolve[vectorEquationStart /. {α -> ξ}, {W1, W2, W3}, {time, 0, 10}, {x, 0, 2 π}, {y, 0, 2 π}, MaxStepFraction -> 1/25]


And then, for looking at the answers:

RW = W /. soln[[1]];

pics = Table[Labeled[Plot3D[RW[[1]] /. {z -> 1}, {x, 0, 2 π}, {y, 0, 2 π}], Row[{"t = ", Pane[time, {50, 12}]}]], {time, 0, 10, .5}]; ListAnimate[pics]

-

## closed as off-topic by Louis, MarcoB, Yves Klett, Mr.Wizard♦Feb 16 at 17:00

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – Louis, MarcoB, Yves Klett, Mr.Wizard
If this question can be reworded to fit the rules in the help center, please edit the question.

It's hard to say without to see the exact problem you try to solve. When you use a StepMonitor or EvaluationMonitor what does that show? Does it get stuck in a specific region? Have you tried to lower AccuracyGloal and PrecisionGoal. If MaxStepFraction is not necessarily the problem leave it away. Are initial/boundary conditions machine numbers? .... – user21 Jun 25 '13 at 7:05
I'm new enough I didn't even know what StepMonitor was. I learned just enough to set this up. You can see exactly what I've typed for initial data/boundary data above. Essentially I just call the three W functions some basic trig function (though it's not important what) and then I use periodic boundary conditions (which is important). The first run finally finished. It was messed up though. Within a second and a half the value for W1 was around 10^50, which is just silly. It's always stabilized around 1 or less before. – James Dilts Jun 26 '13 at 16:15
It's unlikely that you will get an answer; post an example that can be copied to a notebook and run with additional setup. – user21 Jun 26 '13 at 17:31
Could you specify tensor ranks? As I understand, $W$ is $(0,1)$, $(LW)$ is $(0,2)$, and $S$ is $(0,3)$? Also, is the last term $\frac 2n (\mathop{\mathrm{div}} W) g_{ij}$? And your torus is $S^1 \times S^1 \times S^1$? If so, why do you use Cartesian coordinates? Your equation will probably simplify significantly if you write it in terms of angles, no? – Akater May 20 '14 at 5:10