Solve numerically 2D parametric boundary problem

For a personal project, I want to solve the following problem:

Let $$f\left(\boldsymbol{r}\right)$$ be a 2D function, and let $$\boldsymbol{r}\left(s\right)=\left(x\left(s\right),y\left(s\right)\right)$$ be a parameterization of a path in 2D space.
I want to solve the following (horrible) ODE: $$0=\nabla f\left(\boldsymbol{r}\left(s\right)\right)+\frac{\partial\boldsymbol{r}}{\partial s}\frac{\partial f\left(\boldsymbol{r}\left(s\right)\right)}{\partial s}+\left(\frac{\partial^{2}\boldsymbol{r}}{\partial s^{2}}-\left|\frac{\partial^{2}\boldsymbol{r}}{\partial s^{2}}\right|^{2}\frac{\partial\boldsymbol{r}}{\partial s}\right)f\left(\boldsymbol{r}\left(s\right)\right)$$ The solution should be a path $$\boldsymbol{r}\left(s\right)$$, that is defined by a starting point and an ending point $$\boldsymbol{r}\left(0\right)=\boldsymbol{r}_i,\boldsymbol{r}\left(1\right)=\boldsymbol{r}_f$$.

This problem is too hard (if not completely impossible) analytically, so I tried to use Mathematica to solve it numerically for specific cases.

For simple cases such as $$f\left(\boldsymbol{r}\right)=5x$$ Mathematica works fantastic. The problem is that for more complex scenarios, such as $$f\left(\boldsymbol{r}\right)=e^{-\left(x-4\right)^2-y^2}$$, Mathematica gets stuck and never gets to a solution, or returns $$1/0$$ errors.

This is the code I wrote to solve the problem:

(*Definitions*)
ri = {1, 1};
rf = {7, 1};
r[s_] := {x[s], y[s]};
f[r[s]] = Exp[-(x[s] - 4)^2 - y[s]^2];

(*Equations*)
LHS = Grad[f[r[s]], r[s]] + r'[s]*D[f[r[s]], s] + (r''[s] - Norm[r''[s]]^2*r'[s])*f[r[s]];
RHS = {0, 0};
eqs = Table[LHS[[i]] == RHS[[i]], {i, 1, Length[LHS]}]

(*Boundary Conditions*)
conds =
Flatten@Table[{r[0][[i]] == ri[[i]], r[1][[i]] == rf[[i]]}, {i, 1, Length[r[s]]}]

(*Solve and plot*)
sol = NDSolve[Flatten[{eqs, conds}], r[s], {s, 0, 1}];
ParametricPlot[Evaluate[{x[s], y[s]} /. sol], {s, 0, 1}]


My question is: Is there a way to write the code smarter and more efficiently, so that it would work on more complex and non-trivial cases?

Of course, I would love to send any other information you may like to solve the problem.

Thanks a lot!

• Are you looking for a real solution or complex one? Commented Feb 4, 2021 at 12:57
• Real solution @AlexTrounev Commented Feb 4, 2021 at 14:12
• Last solution I got with colocation points 64 shows that there are no unique numerical solution without additional restrictions. Commented Feb 5, 2021 at 18:30

For real solutions we can use collocation method with Haar wavelets. The problem with this method is that we can get one (optimal) solution, but we can't get second if it exists. In the end we can test code on the known (from NDSolve) solution for $$f(r)=5 x$$. First we define equations, solutions, colocation points and wavelets as follows

(*Definitions*)ri = {1, 1};
rf = {7, 1};
r[s_] := {x[s], y[s]};
f[r[s]] = Exp[-(x[s] - 4)^2 - y[s]^2];

(*Equations*)
r'[s]*D[f[r[s]], s] + (r''[s] - (r''[s].r''[s])*r'[s])*f[r[s]];
(*Boundary Conditions*)
conds = Flatten@
Table[{r[0][[i]] == ri[[i]], r[1][[i]] == rf[[i]]}, {i, 1, Length[r[s]]}];

(*Solve and plot*)

eqn = LHS /. {x -> f0, x' -> f1, x'' -> f2, y -> t0, y' -> t1, y'' -> t2}

A = 0; B = 1; J = 3; M = 2^J; dx = (B - A)/(2*M); h1[x_] := Piecewise[{{1, A <= x <= B}, {0, True}}]; p1[x_, n_] := (1/n!)*(x - A)^n;
h[x_, k_, m_] := Piecewise[{{1, Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]}, {-1, Inequality[(1 + 2*k)/(2*m), LessEqual, x, Less, (1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := Piecewise[{{0, x < k/m}, {(-(k/m) + x)^n/n!, Inequality[k/m, LessEqual, x, Less, (1 + 2*k)/(2*m)]}, {((-(k/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!,
(1 + 2*k)/(2*m) <= x <= (1 + k)/m}, {((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 2*(-((1 + 2*k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0]
xl = Table[A + l*dx, {l, 0, 2*M}]; xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2*M + 1}];
f2[x_] := Sum[af[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*h1[x];
f1[x_] := Sum[af[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*p1[x, 1] + f20;
f0[x_] := Sum[af[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0*p1[x, 2] + f20*x + f10;
t2[x_] := Sum[at[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + c0*h1[x];
t1[x_] := Sum[at[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + c0*p1[x, 1] + t10;
t0[x_] := Sum[at[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + c0*p1[x, 2] + t10*x + t00;
bc0 = {f0[0] - 1 == 0, t0[0] - 1 == 0};
bc1 = {f0[1] - 7 == 0, t0[1] - 1 == 0};
var = Flatten[Table[{af[i, j], at[i, j]}, {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
cons = Join[bc0, bc1];

eq =  Flatten[Table[eqn, {s, xcol}]];


Now we optimize solution with NMinimize

varX = Join[{a0, c0, f10, f20, t00, t10}, var]; solM1 =
NMinimize[{Norm[eq], cons}, varX]

(*Out[]= {1.17926*10^-21, {a0 -> 119.195, c0 -> -184.494, f10 -> 1.,
f20 -> -234.87, t00 -> 1., t10 -> 54.8606, af[0, 0] -> 661.782,
at[0, 0] -> 152.644, af[0, 1] -> 166.328, at[0, 1] -> 160.031,
af[1, 1] -> 143.127, at[1, 1] -> -196.665, af[0, 2] -> -109.997,
at[0, 2] -> 177.862, af[1, 2] -> 7.81603, at[1, 2] -> -22.2762,
af[2, 2] -> 6.24531, at[2, 2] -> 12.8123, af[3, 2] -> -58.2248,
at[3, 2] -> -111.313, af[0, 3] -> -76.8545, at[0, 3] -> 232.461,
af[1, 3] -> -20.4284, at[1, 3] -> -154.635, af[2, 3] -> -25.7341,
at[2, 3] -> 167.534, af[3, 3] -> -17.3964, at[3, 3] -> -44.6488,
af[4, 3] -> 9.4494, at[4, 3] -> -3.34721, af[5, 3] -> -106.867,
at[5, 3] -> -50.5085, af[6, 3] -> -4.46456, at[6, 3] -> 143.771,
af[7, 3] -> -40.5735, at[7, 3] -> -131.069}}*)


Actually we have 1.17926*10^-21 for norm of equations in colocation points. Therefore it is numerical solution, and we can use it to plot x, y as

Plot[Evaluate[{f0[s], t0[s]} /. solM1[[2]]], {s, 0, 1},
PlotLegends -> {"x", "y"}, AxesOrigin -> {0, 0}]


Now we can use FindRoot to improve solution as follows

solF1 = FindRoot[Join[Table[eq[[i]] == 0, {i, Length[eq]}], cons],
Table[{varX[[i]], varX[[i]] /. solM1[[2]]}, {i, Length[varX]}],
MaxIterations -> 1000, WorkingPrecision -> 30] // Quiet;

Plot[Evaluate[{f0[s], t0[s]} /. solF1], {s, 0, 1},
PlotLegends -> {"x", "y"}, AxesOrigin -> {0, 0}]


It is not much differ from the plot above, but FindRoot send a message about very small numbers during the evaluation (why NDSolve complains about 1/0 error). Finally we can compare this method solution (red points) with NDSolve solution (solid lines) for f[r[s]] = 5 x[s]

There is second solution of this problem, and we can get it with number of colocation points 64 (J=5). In this case we have 0 for norm of equations in colocation points. In the picture below shown 4 lines for x, y computed with 2 methods.

• It appears that it works! thank you very much :) Commented Feb 5, 2021 at 23:17
• You are welcome! Commented Feb 6, 2021 at 12:21