# How to solve the problem with the error in FindRoot for non-linear equation?

Whoever is reading this have to be patient as the code is slightly long and I have explain each section in order to explain the problem. So I want to solve the set of non-linear algebraic equation using FindRoot. Here is the first section of the code,

np = 11;
zbound[\[Delta]_] := Re[z /. FindRoot[Cos[\[Delta]] Cosh[3 \[Delta]] + Sin[\[Delta]] Sinh[\[Delta]] + (z(-4\[Delta] Cos[\[Delta]] Cosh[3 \[Delta]] + Sin[\[Delta]] (Cosh[\[Delta]] - 4 \[Delta] Sinh[\[Delta]]) + Cos[\[Delta]] Sinh[3 \[Delta]]))/\[Delta] == 0, {z, 0.1}]]
zvalue = Table[z[i], {i, 1, np}];
guess = Table[(-((np + 1)/2) + i)*(0.5*I)/10000, {i, 1, np}]
eqn[\[Delta]_] := Table[(E^(-I*z[i]*20000) == \!$$\*UnderoverscriptBox[\(\[Product]$$, $$j = 1$$, $$np$$]$$If[j != i, \*FractionBox[\((I*\((\*FractionBox[\(z[i] - z[j]$$, $$2$$])\) - zbound[\[Delta]])\), $$(\(-I$$*$$(\*FractionBox[\(z[i] - z[j]$$, $$2$$])\) - zbound[\[Delta]])\)], 1]\)\)), {i, np}];


Just copy and paste the above code into Mathematica Notebook. The last line is the non-linear equation I want to solve.

eqnsol[\[Delta]_, initialvalue_] :=zvalue /.FindRoot[eqn[\[Delta]],Transpose[{zvalue, initialvalue}], MaxIterations -> 400]


This will find the solution of the equations that depends on $$\delta$$. With the initial value as "guess",

Block[{\[Delta] = 1.60919303}, eqnsol[\[Delta], guess]]


This will give me the solution of a particular value of $$\delta$$. Now, I want to evolve the solution by increasing the value of $$\delta$$.

\[Delta]list = SetPrecision[Range[1.60919303, 1.60919305, 10^-9], 10];
approxData = {};
Do[sol = eqnsol[\[Delta], guess];
AppendTo[approxData, sol];
guess = sol;, {\[Delta], \[Delta]list}]


This will evolve the solution by updating the guess as previous solution. I can plot the data with respect to $$\delta$$.

GraphicsRow[{ListPlot[Table[Transpose[{\[Delta]list, Re[approxData[[All, i]]]}], {i, 1, np}]],ListPlot[Table[Transpose[{\[Delta]list, Im[approxData[[All, i]]]}], {i, 1, np}]]}] No problem till this point. Now if I want to increase the range of $$\delta$$,

\[Delta]list = SetPrecision[Range[1.60919303, 1.6091931, 10^-9], 10];


Running the same loop, you get error, and it spits out some junk values.  Reasoning behind it I believe is that as the solution increases, the exponent in the LHS of the equation is decreasing, and at some point it becomes smaller that the precision of the code. I tried increasing the WorkingPrecision but it didn't work.

Note: If someone is able to solve the problem, try also increasing the value of np to some odd integer to see the solution holds or not.

• @Bill I tried that, but I am not seeing any syntax error. Mar 23 at 21:19
• It is not necessary to use a numeric technique (i.e., FindRoot) in the definition of zbound. Use zbound[\[Delta]_] = Re[SolveValues[Cos[\[Delta]] Cosh[3 \[Delta]] + Sin[\[Delta]] Sinh[\[Delta]] + (z (-4 \[Delta] Cos[\[Delta]] Cosh[3 \[Delta]] + Sin[\[Delta]] (Cosh[\[Delta]] - 4 \[Delta] Sinh[\[Delta]]) + Cos[\[Delta]] Sinh[3 \[Delta]]))/\[Delta] == 0, z][]]; Mar 23 at 22:20

As Bob stated, can solve for $$z$$ explicitly in zbound (I changed the variable names to distinguish my code):

newZBound[\[Delta]_] =
z /. Solve[
Cos[\[Delta]] Cosh[3 \[Delta]] +
Sin[\[Delta]] Sinh[\[Delta]] + (z (-4 \[Delta] Cos[\[Delta]] \
Cosh[3 \[Delta]] +
Sin[\[Delta]] (Cosh[\[Delta]] -
4 \[Delta] Sinh[\[Delta]]) +
Cos[\[Delta]] Sinh[3 \[Delta]]))/\[Delta] == 0, z][[1, 1]]


Now try $$np=25$$:

    np = 25;
newEqn[\[Delta]_] := Table[(E^(-I*z[i]*20000) == \!$$\*UnderoverscriptBox[\(\[Product]$$, $$j = 1$$, $$np$$]$$If[j != i, \*FractionBox[\((I*\(( \*FractionBox[\(z[i] - z[j]$$, $$2$$])\) -
newZBound[\[Delta]])\), $$(\(-I$$*$$( \*FractionBox[\(z[i] - z[j]$$, $$2$$])\) - newZBound[\[Delta]])\)],
1]\)\)), {i, np}];


Now, for the guess values you initially used a pure imaginary number and specified 0.5 which immediately drops the precision to machine precision. Running that by itself, when the number of equation increase above 11, FindRoot encounters a singular Jacobian and recommends perturbing the intitial values. So I first change 0.5 to 1/2 to force arbitrary precision and add a perturb value of 1/1000 in the guess list:

perturbVal = 1/1000;
zvalue = Table[z[i], {i, 1, np}];
guess = Table[(-((np + 1)/2) + i)*(1/2*I)/10000 + perturbVal, {i, 1,
np}]

newEqnSol[\[Delta]_, initialvalue_] :=
zvalue /.
FindRoot[newEqn[\[Delta]], Transpose[{zvalue, initialvalue}],
MaxIterations -> 400]


Now increase the range of $$\delta$$ and plot the resutls:

\[Delta]list = SetPrecision[Range[1.60919303, 1.6091931, 10^-9], 10]
approxData = {};
Do[sol = newEqnSol[\[Delta], guess];
AppendTo[approxData, sol];
guess = sol;, {\[Delta], \[Delta]list}]
GraphicsRow[{ListPlot[
Table[Transpose[{\[Delta]list, Re[approxData[[All, i]]]}], {i, 1,
np}], PlotStyle -> PointSize[0.0075]],
ListPlot[
Table[Transpose[{\[Delta]list, Im[approxData[[All, i]]]}], {i, 1,
np}], PlotStyle -> PointSize[0.0075]]}] • What version do you run? Mar 28 at 16:04
• @josh Awesome solution. However, I don't think this is correct solutions. Yes, the numerical code is running, but if you notice, the sum of all the z's have to be zero. That's why I took purely imaginary guess. This can be checked by taking np=3 and just multiply three equations. You will get $e^{i(z_1+z_2+z_3)*(20000)}=1$. Anyway, this is a good starting point for me to change the guess values and see if I get near my solutions. Mar 28 at 17:32
• @Alex Trounev: I'm running 13.2
– josh
Mar 28 at 17:38