The most common method of creating patterns on a membrane involves drawing geodesic lines on the membrane. Justification of the layout of the geodesic lines is beyond the scope of this text .

Although the geodesic lines can be drawn freely to form patterns of any shape and position, it is common that the route of the geodesic lines follow some models determined:

how to draw these roughly parallel geodesic lines is something that every user usually does in a particular way depending on different objectives I need repeat the same exercise with a geodesic line located on the Plot3D of f, like generatrix line.

For example, let the function

Plot3D[x + y - x y, {x, -5, 5}, {y, -5, 5}, PlotTheme -> "Scientific"]

enter image description here Can I draw their lines of geodesic as shown below

A method of multi-geodesic uses end points of the new geodesic lines found at the edge of the membrane. If we are with a membrane formed by different parts, we can not perform the function of multi-geodesic directly. enter image description here

What we must do is break down the membrane in two simple membranes and find the pattern in each. I would like to create images based on above method as shown below: enter image description here enter image description here Is this possible in Mathematica? I do not know how this can be implemented. Please Help.

  • 1
    $\begingroup$ I don't think "geodesic" means what you think it means here. Anyway, you mentioned "the route of the geodesic lines follow some models determined"; there should be a mathematical description of said "model", no? $\endgroup$ Aug 26, 2017 at 23:05
  • $\begingroup$ What do you mean by geodesic? I don't get what you aim for... $\endgroup$ Aug 27, 2017 at 0:24
  • $\begingroup$ @HenrikSchumacher see this wintess.com/multi-geodesicas-1 $\endgroup$ Aug 27, 2017 at 0:26

1 Answer 1


Plot the surface.

surface = Plot3D[x + y - x y, {x, -5, 5}, {y, -5, 5}, AxesLabel -> {"x", "y", "z"}]


The aim is to draw geodesics - i.e. minimum distance paths - on this surface.

Define the geodesic as a parametric curve.

geo[t_] = {x[t], y[t], z[t]};

Define z[t] in terms of x[t] and y[t].

z[t_] = x[t] + y[t] - x[t] y[t];

A minimum distance path can be computed as the solution of the variational problem of finding a path on the surface whose length is stationary with respect to perturbations of the path.

Get the variational methods package. At this point it is useful to review the Variational Methods documentation.

<< VariationalMethods`

Construct the function that computes the rate of movement along the geodesic.

f[t] = Sqrt[#.#] &[D[{x[t], y[t], z[t]}, t]]

Construct the Euler equations which the geodesic solution must satisfy.

diffeqns = EulerEquations[f[t], {x[t], y[t]}, t]

<omitted messy Euler equations>

Define some initial conditions at t=0 to fix the position and direction at the start of the geodesic.

initconds = {x[0] == -5, y[0] == 5, x'[0] == 1, y'[0] == -1};

Define the maximum value of the geodesic parameter to be used. This is used to determine how far to compute the numerical solution of the Euler equations.

tmax = 9;

Numerically solve the Euler equations subject to the initial conditions.

This throws the warning "NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.", but this appears to be a benign problem.

EDIT: I wrapped diffeqns in FullSimplify to fix an error that occurred in version 12.0.

solution = NDSolve[Join[FullSimplify[diffeqns], initconds], {x, y}, {t, 0, tmax}];

Plot the geodesic, and superimpose it on the plot of the surface.

geodesic = ParametricPlot3D[
  {x[t], y[t], z[t]} /. solution[[1]] // Evaluate, {t, 0, tmax}];
Show[surface, geodesic]

enter image description here

This geodesic is parameterised by the initial conditions and tmax, and these parameters can be varied to generate other geodesic solutions.

WARNING: I found that the numerical solution of the Euler equations is not in general as easy to obtain as in the specific example above, which is therefore merely a guide to how variational methods can be used to compute geodesics subject to initial conditions. More generally, you will have to tweak the NDSolve options, and so on.

  • $\begingroup$ But running your code in mma12.0 will return an error. $\endgroup$ Feb 9, 2020 at 9:54
  • 2
    $\begingroup$ I fixed the problem in NDSolve by wrapping diffeqns in FullSimplify. $\endgroup$ Apr 9, 2020 at 21:58
  • $\begingroup$ @stephen, would you please provide the complete code? $\endgroup$ Jul 3, 2022 at 18:03
  • $\begingroup$ I just copied (to Mathematica 13.1) and ran each line of code listed in my answer above, and everything works as described. The "WARNING" (at the end of my answer) about tweaking options was based on some "scratch" computations that have no formal record in my notebook archive. $\endgroup$ Jul 3, 2022 at 18:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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