This problem can be solved with collocation method using Bernoulli wavelets. We don't need to run NDSolve
and NMinimize
so it is alternative to that proposed by Sjoerd Smit. Nevertheless we can compare final result for npts=8
. First we define wavelets, functions and derivatives as follows
n = 2;
M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1; U = Array[u, {M}];
dx = 1/M; A = 0; xl = Table[A + l*dx, {l, 0, M}]; tcol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}];
psi1[x_] := Piecewise[{{BernoulliB[2, x], 0 <= x < 1}, {0, True}}];
psi2[x_] := Piecewise[{{BernoulliB[1, x], 0 <= x < 1}, {0, True}}];
psi1jk[x_, j_, k_] := psi1[j*x - k];
psi2jk[x_, j_, k_] := psi2[j*x - k];
psijk = Compile[{{x, _Real}, {j, _Integer}, {k, _Integer}}, (psi1jk[x,
j, k] + psi2jk[x, j, k])/2];
(*intjk=Integrate[psijk[x,j,k],x,Assumptions\[Rule]{j>0,k>0}]*)
psijk1 = Compile[{{x, _Real}, {j, _Integer}, {k, _Integer}},
Piecewise[{{(-k + k^2)/(2*j),
j > 0 && k == 0 &&
1/j - x < 0}, {(1/6)*(-x + 3*k^2*x - 3*j*k*x^2 + j^2*x^3),
j > 0 && k == 0 && x > 0 &&
1/j - x >=
0}, {(k - k^3 - j*x + 3*j*k^2*x - 3*j^2*k*x^2 + j^3*x^3)/(6*
j), j > 0 && k > 0 && k/j - x < 0 && 1/j + k/j - x >= 0}},
0]];
Psi[x_] :=
Join[{1},
Flatten[Table[psijk[x, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
int1[x_] :=
Join[{x},
Flatten[Table[
psijk1[x, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
var1 = Join[{a0},
Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var2 =
Join[{b0},
Flatten[Table[b[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var3 =
Join[{c0}, Flatten[Table[c[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
z1[t_] := var3.Psi[t]; z[t_] := var3.int1[t] + c1;
y1[t_] := var1.Psi[t]; y[t_] := var1.int1[t] + a1;
x1[t_] := var2.Psi[t]; x[t_] := var2.int1[t] + b1;
Second, we define model to be optimized
varM = Join[{a1, b1, c1}, var1, var2, var3]; tf = 1;
j = 0;
AbsoluteTiming[Do[j = j + 1; eq = Flatten[Table[{2*x[xcol[[i]]] - 3*u[i]*x[xcol[[i]]] - x1[xcol[[i]]]/tf - x[xcol[[i]]]*y[xcol[[i]]] == 0,
(-u[i])*x[xcol[[i]]] + y[xcol[[i]]] - 2*x[xcol[[i]]]*y[xcol[[i]]] - y1[xcol[[i]]]/tf == 0,
-z1[xcol[[i]]]/tf + (x[xcol[[i]]] - 1)^2 + (y[xcol[[i]]] - 1)^2 == 0}, {i, Length[xcol]}]];
cons = Join[eq, {x[0] == 1, y[0] == 1/2, z[0] == 0}]; sol[j] = Quiet[FindRoot[cons, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]];
s[j] = z[1] /. sol[j]; uu[j] = U; , {u[1], 0, 1, 1}, {u[2], 0, 1, 1}, {u[3], 0, 1, 1}, {u[4], 0, 1, 1}, {u[5], 0, 1, 1}, {u[6], 0, 1, 1},
{u[7], 0, 1, 1}, {u[8], 0, 1, 1}]]
Third, we can fined minimum and plot solution as follows
{km, sm} = MinimalBy[Table[{k, s[k]}, {k, 1, j}], Last][[1]]
(*{227, 0.786878}*)
lst1 = Join[{{0, 1}},
Table[{xcol[[i]], x[xcol[[i]]] /. sol[km]}, {i, M}]]; lst2 =
Join[{{0, 1/2}}, Table[{xcol[[i]], y[xcol[[i]]] /. sol[km]}, {i, M}]];
ListLinePlot[{lst1, lst2}, PlotLabel -> uu[km], Frame -> True,
Axes -> False, PlotLegends -> {"x", "y"}]

So far so good, we have practically same minimum value 0.786878
compare to solution by Sjoerd Smit of 0.787574
.
Now we put n=3
in my code that corresponds to npts = 16
. In this case model to be optimized is
varM = Join[{a1, b1, c1}, var1, var2, var3]; tf = 1;
j = 0;
AbsoluteTiming[Do[j = j + 1; eq = Flatten[Table[{2*x[xcol[[i]]] - 3*u[i]*x[xcol[[i]]] - x1[xcol[[i]]]/tf - x[xcol[[i]]]*y[xcol[[i]]] == 0,
(-u[i])*x[xcol[[i]]] + y[xcol[[i]]] - 2*x[xcol[[i]]]*y[xcol[[i]]] - y1[xcol[[i]]]/tf == 0,
-z1[xcol[[i]]]/tf + (x[xcol[[i]]] - 1)^2 + (y[xcol[[i]]] - 1)^2 == 0}, {i, Length[xcol]}]];
cons = Join[eq, {x[0] == 1, y[0] == 1/2, z[0] == 0}]; sol[j] = Quiet[FindRoot[cons, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]];
s[j] = z[1] /. sol[j]; uu[j] = U; , {u[1], 0, 1, 1}, {u[2], 0, 1, 1}, {u[3], 0, 1, 1}, {u[4], 0, 1, 1}, {u[5], 0, 1, 1}, {u[6], 0, 1, 1},
{u[7], 0, 1, 1}, {u[8], 0, 1, 1}, {u[9], 0, 1, 1}, {u[10], 0, 1, 1}, {u[11], 0, 1, 1}, {u[12], 0, 1, 1}, {u[13], 0, 1, 1},
{u[14], 0, 1, 1}, {u[15], 0, 1, 1}, {u[16], 0, 1, 1}]]
With my code we got final result in 26 min at j=65536
. New result is
{km, sm} = MinimalBy[Table[{k, s[k]}, {k, 1, j}], Last][[1]]
(*Out[]= {63763, 0.782377} *)
This result not so differ from above while picture not looks similar

This code can be optimize with Compile
and ParallelDo
as it shown here. We also can compare last picture with NMimimize
(it finished after 8 h 27 min running on my laptop). Final result looks like
min = NMinimize[{Indexed[sol @@ uvec, 1],
uvec \[Element] Integers && And @@ Map[# == 0 || # == 1 &, uvec]},
uvec] // AbsoluteTiming
Out[]= {30454.7, {0.783045, {u[1] -> 1, u[2] -> 1, u[3] -> 1,
u[4] -> 1, u[5] -> 1, u[6] -> 0, u[7] -> 0, u[8] -> 0, u[9] -> 1,
u[10] -> 0, u[11] -> 0, u[12] -> 1, u[13] -> 0, u[14] -> 0,
u[15] -> 1, u[16] -> 0}}}
Now we have difference in vector U
as well, while minimum value with my code and with NMinimize
- 0.782377 and 0.783045 consequently. Pictures also are differ but it can be explained by different interpolation method with wavelets and NDSolve
.
