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

Question

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.

Answer 1

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]]}]