Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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]
share|improve this 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

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.