# Problem of step size effectively zero

I've been trying to solve the next system of differential equations which is very similar to this one in which I also sought help

Step size is effectively zero

$$F^2-G^2+HF'-F''+1=0$$

$$2GF+HG'-G''=0$$

$$2F+H'=0$$

$$HH'+P'-H''=0$$

With boundary conditions

$$F(0)=G(0)=H(0)=P(0)=0$$

$$F(∞)=0, G(∞)=1$$

Again, I turned them into a system of first order ordinary differential equations, $$F=x, F'=y, G=z, G'=s, H'=p$$ and $$P=u$$. So I tried the same code that run with the original system, with a few corrections, but it doesn't seem to work with this system, probably because it is unstable too. At t = 3.45 there is a step size effectively zero problem. I also tried with ParametricNDSolve, but I really don't know how it works. In the original system I had an idea of the behavior of the solution unlike this case in which I don't, so I assumed P(14) = 0, in an attempt to copy the idea of getting a similar solution and then improve it.

odes = {x'[t] == y[t], y'[t] == x[t]^2 - z[t]^2 + s[t] y[t] + 1,
s'[t] == 2 z[t] x[t] + p[t] s[t], z'[t] == s[t], p'[t] == -2 x[t],
u'[t] == 2 x[t] p[t] - 2 y[t]};
bcs = {x[0] == 0, z[0] == 0, p[0] == 0, x[14] == 0, z[14] == 1,
u[0] == 0};
vars = {x, y, z, s, p, u};
sol32 = NDSolve[{odes, bcs}, vars, {t, 0, 14},
Method -> {"Shooting",
"StartingInitialConditions" -> {s[14] == 0, y[14] == 0,
p[14] == 0, x[14] == 0, z[14] == 1, u[14] == 0}},
PrecisionGoal -> 6, AccuracyGoal -> 6, WorkingPrecision -> 32];

solMP = NDSolve[{odes, bcs}, vars, {t, 0, 14},
Method -> {"Shooting",
"StartingInitialConditions" -> {Through[vars[14]] ==
Through[vars[14] /. First[sol32]]}}, PrecisionGoal -> 10,
AccuracyGoal -> 10]

ListLinePlot[{x, z, p} /. First[sol], PlotLegends -> {F, G, H}]



It is the exact same one in the original post but with the equation $$u'= 2 x p - 2y$$, and $$u(14) = 0$$.

• Are you looking solution with NDSolve only or with any other methods? Commented Mar 2, 2021 at 17:22
• Any method that works will be fine, it is highly unlikely that this system has an analytical solution though. Commented Mar 2, 2021 at 22:02

## 1 Answer

This problem can be solved with collocation method using Bernoulli wavelets. First we map interval to the unit interval and define wavelets, functions and derivatives for vars = {x, y, z, s, p, u} on the unit domain as follows

Clear["Global*"]

n = 4;
M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1;
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}]]]; var4 =
Join[{d0},
Flatten[Table[d[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var5 =
Join[{e0},
Flatten[Table[e[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]]; var6 =
Join[{f0}, Flatten[Table[f[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;
s1[t_] := var4 . Psi[t]; s[t_] := var4 . int1[t] + d1;
p1[t_] := var5 . Psi[t]; p[t_] := var5 . int1[t] + e1;
u1[t_] := var6 . Psi[t]; u[t_] := var6 . int1[t] + f1;


Second, we define boundary conditions and equations to be optimize

L = 14;
bcs = {x[0] == 0, z[0] == 0, p[0] == 0, x[1] == 0, z[1] == 1,
u[0] == 0};

eq = Flatten[Table[{-x1[xcol[[j]]]/L + y[xcol[[j]]], 1 + x[xcol[[j]]]^2 + s[xcol[[j]]]*y[xcol[[j]]] - y1[xcol[[j]]]/L - z[xcol[[j]]]^2,
p[xcol[[j]]]*s[xcol[[j]]] - s1[xcol[[j]]]/L + 2*x[xcol[[j]]]*z[xcol[[j]]], s[xcol[[j]]] - z1[xcol[[j]]]/L,
-p1[xcol[[j]]]/L - 2*x[xcol[[j]]], -u1[xcol[[j]]]/L + 2*p[xcol[[j]]]*x[xcol[[j]]] - 2*y[xcol[[j]]]}, {j, M}]];
varM = Join[{a1, b1, c1, d1, e1, f1}, var1, var2, var3, var4, var5, var6];


Finally we use NMinimize to solve this problem

sol = NMinimize[{Norm[eq], bcs}, varM]


With this bcs we have norm of equations of 3.83121*10^-6 and it is not bad. For visualization we use Plot so Compile complains, it is why we use Plot with Quite

Plot[Evaluate[{x[t/L], z[t/L], p[t/L]} /. sol[[2]]], {t, 0, 14},
PlotLegends -> {"F", "G", "H"}, Frame -> True] // Quiet


We can improve optimal numerical solution using FindRoot as follows

sol1 = FindRoot[Join[Table[eq[[i]] == 0, {i, Length[eq]}], bcs],
Table[{varM[[i]], varM[[i]] /. sol[[2]]}, {i, Length[varM]}],
MaxIterations -> 1000]


But it produces same plot

Plot[Evaluate[{x[t/L], z[t/L], p[t/L]} /. sol1], {t, 0, 14},
PlotLegends -> {"F", "G", "H"}, Frame -> True] // Quiet


For some implementation of the colocation method with Bernoulli wavelets see my post on this forum

High precision numerical solution of the nonlinear Volterra integral equation

There are also answer here with application Bernoulli wavelets to optimization problem. Some papers about this method

A new numerical method for fractional order Volterra integro-differential equations

The Bernoulli wavelets operational matrix of integration and its applications for the solution of linear and nonlinear problems in calculus of variations

• Thanks! I would not have been able to come up with that code, pretty interesting that the solution has remained virtually the same. I run the code and plotted P, what I was looking for, too. Commented Mar 3, 2021 at 16:18
• @SebastiánFrades You are welcome! This is more effective method for BVP then "Shooting" method implemented in NDSolve`. Commented Mar 3, 2021 at 17:05
• Sorry to bother you, I'd really appreciate if could you tell me where to read about collocation method using Bernoulli wavelets, or how to use this method in mathematica because I haven't been able to find much on the internet. Commented Mar 6, 2021 at 23:09
• @SebastiánFrades please, see update to my answer. Commented Mar 6, 2021 at 23:45