1
$\begingroup$

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

enter image description here

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.

enter image description here

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ @Bill I tried that, but I am not seeing any syntax error. $\endgroup$
    – sslucifer
    Mar 23, 2023 at 21:19
  • 1
    $\begingroup$ 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][[1]]]; $\endgroup$
    – Bob Hanlon
    Mar 23, 2023 at 22:20

1 Answer 1

3
+100
$\begingroup$

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

enter image description here

$\endgroup$
3
  • $\begingroup$ What version do you run? $\endgroup$ Mar 28, 2023 at 16:04
  • $\begingroup$ @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. $\endgroup$
    – sslucifer
    Mar 28, 2023 at 17:32
  • $\begingroup$ @Alex Trounev: I'm running 13.2 $\endgroup$
    – josh
    Mar 28, 2023 at 17:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.