# Wrong Boundary Conditions Result Using Wavelet Collocation

I have a functional $$S$$,

$$S = \int_{x_0}^{x_b} dx \frac{1}{z(x)^d} \sqrt{1 + \frac{z'(x)^2}{f(z)}}, \qquad f(z) = 1-\left(\frac{z(x)}{z_h}\right)^{d+1}$$

where $$d=3$$ is the dimension and $$z_h$$ is some constant. We can solve the equation of motion through the Euler-Lagrange equation which is a nonlinear second-order differential equation in this case. The boundary conditions are given by,

$$z(x_0) = z_s = 9.3, \; z(x_b) = \epsilon = 10^{-5}$$

The ODE is solved using wavelet collocation, specifically Haar wavelets. However, after getting the solution and checking if the solution satisfies the boundary points, I found that it did not match the imposed boundary conditions. On the other hand, I also observed that when the factor $$1/z(x)^d$$ is removed from the functional $$S$$, the solution correctly satisfies the boundary condition imposed. See the code below for the situation.

Just to give more information, the reason the wavelet approach was used is because I cannot produce the solution using the shooting method (so many issues, ugh!). Now, at least a solution can be obtained, but it does not match the boundary conditions. Is the factor $$1/z(x)^d$$ some specific problem that should be addressed by some specific technique? Any clue as to why is it like that and how to resolve this?

ClearAll["Global*"]
Needs["VariationalMethods"]
f = 1 - (z[x]/zh)^(d + 1);
L = (Sqrt[1 + (z'[x]^2/f)]/z[x]^d);(*Lagrangian*)
eulageq = EulerEquations[L, z[x], x];(*Euler-Lagrange equation*)
s = Solve[eulageq, z''[x]][[1]] // Simplify;
eq0 = z''[x] == s[[1, 2]];

(*Haar Wavelet*)
J = 3;
M = 2^J;
dx = 1/(2 M);
xl = Table[l*dx, {l, 0, 2 M}];
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2 M + 1}];
h1[x_] := Piecewise[{{1, 0 <= x <= 1}, {0, True}}];
p1[x_, n_] := (1/n!) x^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]
var1 = Flatten[Table[a[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
z2[x_] := Sum[a[i, j] h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0 h1[x];
z1[x_] := Sum[a[i, j] p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0 p1[x, 1] + a1;
z0[x_] := Sum[a[i, j] p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + a0 p1[x, 2] + a1 x + a2;
eqs = {eq0} /. {z''[x] -> z2[x], z'[x] -> z1[x], z[x] -> z0[x]};
var = Join[var1, {a0, a1, a2}];
x0 = 10^-10;
xb = 10^-1;
y0 = (x0 - x0)/(xb - x0);
yb = (xb - x0)/(xb - x0);
eq = Table[eqs /. {d -> 3, zh -> 10}, {x, xcol}] // Flatten;
bc = {z0[y0] == 93/10, z0[yb] == 10^-5};

sol = FindRoot[Join[eq, bc], Table[{var[[i]], 9/10}, {i, Length[var]}], Method -> {"Newton", "StepControl" -> "TrustRegion"}, MaxIterations -> 10^6, WorkingPrecision -> 20];

Evaluate[z0[y0] /. sol]
8.0495553936380335557

Evaluate[z0[yb] /. sol]
2.883219672021461112

(*1/z[x]^d removed*)
Evaluate[z0[y0] /. sol]
9.3000000000000000000

Evaluate[z0[yb] /. sol]
0.000010000000000000

The code is written in Mathematica. Credits to Alex Trounev (member of Mathematica Stack) who wrote the Haar wavelet code.

• It could be better to show equation eq0 in Latex form as follows $z''(x)-\frac{-2 d \text{zh}^2 \left(z'(x)^2+1\right)+\text{zh} z(x) \left(\frac{z(x)}{\text{zh}}\right)^d \left((d-1) z'(x)^2+4 d\right)-2 d z(x)^2 \left(\frac{z(x)}{\text{zh}}\right)^{2 d}}{2 \text{zh} z(x) \left(\text{zh}-z(x) \left(\frac{z(x)}{\text{zh}}\right)^d\right)}=0$ Commented Dec 9, 2022 at 15:54

This problem can be solved with the Shooting method for inverse function $$x(z)$$ as follows

ClearAll["Global*"]
Needs["VariationalMethods"]
f = 1 - (z/zh)^(d + 1);
L = (Sqrt[1 + (x'[z]^(-2)/f)]/z^d) x'[z];(*Lagrangian*)
eulageq = EulerEquations[L, x[z], z];(*Euler-Lagrange equation*)
s = Solve[eulageq, x''[z]][[1]] // Simplify;
eq0 = x''[z] - s[[1, 2]] == 0

Please, note, that eq0 has a form

$$\frac{x'(z) \left(2 d \left(\text{zh}-z \left(\frac{z}{\text{zh}}\right)^d\right)^2 x'(z)^2+\text{zh} \left(2 d \text{zh}-(d-1) z \left(\frac{z}{\text{zh}}\right)^d\right)\right)}{2 z \text{zh} \left(z \left(\frac{z}{\text{zh}}\right)^d-\text{zh}\right)}+x''(z)=0$$

BVP solution

eq1 = eq0 /. {zh -> 10, d -> 3};
zs = 93/10; z0 = 10^-5; sol =
NDSolveValue[{eq1, x[zs] == 10^-10, x[z0] == 1/10}, x, z]

Visualization $$x(z)$$ and $$z(x)$$

Plot[sol[z], {z, z0, zs}, AxesLabel -> {"z", "x"}, PlotRange -> All]
lst = Table[{sol[z], z}, {z, z0, zs, .01}]; ListLinePlot[lst,
InterpolationOrder -> 2, Frame -> True, FrameLabel -> {"x", "z"}]

We also can solve this problem by using Haar wavelets as follows

J = 5; M = 2^J; dx = 1/(2*M); xl = Table[l*dx, {l, 0, 2*M}];
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2*M + 1}];
h1[x_] := Piecewise[{{1, 0 <= x <= 1}, {0, True}}];
p1[x_, n_] := (1/n!)*x^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]
var1 = Flatten[Table[a[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
X2[x_] :=
Sum[a[i, j]*h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0*h1[x];
X1[x_] :=
Sum[a[i, j]*p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0*p1[x, 1] + a1;
X0[x_] :=
Sum[a[i, j]*p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] +
a0*p1[x, 2] + a1*x + a2;

xb = 10^-1; zb = zs; z0 =
10^-5; eq2 = {eq1} /. {x''[z] -> X2[Z]/zb^2, x'[z] -> X1[Z]/zb,
z -> zb Z};
var = Join[{a0, a1, a2}, var1];
x0 = 10^-10;
y0 = z0/zs;
yb = 1;
eq = Table[eq2, {Z, xcol}] // Flatten;
bc = {X0[y0] == xb, X0[yb] == x0};

sol1 = FindRoot[Join[eq, bc],
Table[{var[[i]], 1/10}, {i, Length[var]}]];

Visualization sol1 (red line) together with sol (blue line)

The problem arises when we try to solve the equation for $$z(x)$$. The corresponding equation shown in my comment and for $$d=3, zh=10$$ is given by

$$z''(x)-\frac{\frac{1}{100} z(x)^4 \left(2 z'(x)^2+12\right)-600 \left(z'(x)^2+1\right)-\frac{3 z(x)^8}{500000}}{20 z(x) \left(10-\frac{z(x)^4}{1000}\right)}=0$$

We use the Haar wavelets method described above and adapted for this case as follows

eqZ = Derivative[2][z][x] - (-(2*d*zh^2*(Derivative[1][z][x]^2 + 1)) +

zh*z[x]*(z[x]/zh)^d*((d - 1)*Derivative[1][z][x]^2 + 4*d) -
2*d*z[x]^2*(z[x]/zh)^(2*d))/(2*zh*
z[x]*(zh - z[x]*(z[x]/zh)^d)) /. {d -> 3, zh -> 10}

eqZ2 = eqZ /. {z''[x] -> X2[Z]/xb^2, z'[x] -> X1[Z]/xb, z[x] -> X0[Z]};

bc = {X0[1] == z0, X0[x0/xb] == zs};
solz = FindRoot[
Join[Table[eqZ2 == 0 /. Z -> xcol[[i]], {i, Length[xcol]}] //
Flatten, bc], Table[{var[[i]], 1/10}, {i, Length[var]}],
MaxIterations -> 10000, WorkingPrecision -> 30,
Method -> {"Newton", "StepControl" -> "TrustRegion"}];

If we visualize the NDSolve solution (blue line), the expected solution (red points), and solz (black line) then we see large discrepancies in the two solutions.

Note, that the same solution we compute also when the initial guess is given by red points. Obviously, this error is related to the method used to find the roots. We can discuss this issue with a company that produces Mathematica. But beyond that, how can we solve this problem?