3
$\begingroup$

Suppose we have set of matrix equations, for this case lets consider a 2x2 set. The equations are:

α*x*840.92 == 470.44*(Coth[α] - 1/α)+1329.74(Coth[β] - 1/β)

β*x*840.92 == 664.87*(Coth[α] - 1/α)+143.61(Coth[β] - 1/β)

I have tried using NSolve to no avail as it keeps the Coth unevaluated.

x=0.1
NSolve[{α*x*840.92 == 470.44*(Coth[α] - 1/α) + 1329.74*(Coth[β] - 1/β),β*x *840.92 == 664.87*(Coth[α] - 1/α) + 143.61*(Coth[β] - 1/β) && α > 0 && β > 
0}, {α, β}, Reals]

This just returns the same message.

The main idea is that for a given x, α and β can be calculated. I was trying to use a series expansion but the results are very unreliable I think. This should be able to work for any size matrix.

| improve this question | | | | |
$\endgroup$
  • 1
    $\begingroup$ For those contemplating (0,0) as a solution, note the system is undefined at (0,0). You'd have to clear denominators (including Coth) first, and then you can just plug in (0,0) to set it's a solution of the resulting system. $\endgroup$ – Michael E2 Jan 11 '18 at 2:14
7
$\begingroup$

Updated for clarity

It turns out there's a branch of solutions for 0 < x <= 0.500003 that you can get at with FindRoot. First, define

eq1 = α*x*840.92 == 470.44*(Coth[α] - 1/α) + 1329.74 (Coth[β] - 1/β);
eq2 = β*x*840.92 == 664.87*(Coth[α] - 1/α) + 143.61 (Coth[β] - 1/β);

and then rearrange these equations to get two definitions for x as a function of α and β:

x1[α_, β_] := (470.44*(Coth[α] - 1/α) + 1329.74 (Coth[β] - 1/β)) / (α*840.92)
x2[α_, β_] := (664.87*(Coth[α] - 1/α) + 143.61 (Coth[β] - 1/β)) / (β*840.92)

Then for any given x, the triplet {α, β, x} solves eq1 && eq2 if and only if x1[α, β] == x2[α, β] == x.

Finding the solutions

To get a sense of what's going on, we can plot x1 and x2 over the {α, β}-plane:

Plot3D[{x1[α, β], x2[α, β]}, {α, 0, 2}, {β, 0, 1}, AxesLabel -> {"α", "β", "x"}, 
 PlotRange -> {All, All, {0, 1}}]

enter image description here

So we're looking for the intersection between those two surfaces. It's hard to see in the figure without rotating it, but there's something weird going on around x == 0.5. I'll go into a little more detail on this below, but basically, if you look at it as a bifurcation problem with x as the bifurcation parameter, then it turns out there's a nontrivial branch of solutions for x <= 0.500003.

Getting back to a numerical search for the solutions, you can use

x0 = 0.500003;
αβsol[x_ /; 0 < x <= x0] := FindRoot[
  {x1[α, β] == x, x2[α, β] == x}, {{α, 2., 10^-10, ∞}, {β, 1., 10^-10, ∞}}
]

or, if you just want the values and not replacement rules

αβsolvalue[x_ /; 0 < x <= x0] := {α, β} /. FindRoot[
  {x1[α, β] == x, x2[α, β] == x}, {{α, 2., 10^-10, ∞}, {β, 1., 10^-10, ∞}}
]

(where I've restriced α and β to be >= 10^-10, and started every search at {α, β} == {2, 1} simply because it works).

Plotting the solutions:

Plot[{αβsolvalue[x][[1]], αβsolvalue[x][[2]]}, {x, 0, x0}, 
  PlotLegends -> {α, β}, AxesLabel -> {x, None}]

enter image description here

Note that {α, β} -> 0 as x -> x0 and that the solution is undefined for x = 0.

Adding these solutions to the above plot of x1[α, β] and x2[α, β] shows that it's found the intersection of those two planes:

Show[
 Plot3D[
   {x1[α, β], x2[α, β]}, {α, 0, 2}, {β, 0, 1}, 
   PlotRange -> {{0, 2}, {0, 1}, {0, 1}}],
 ParametricPlot3D[
   Flatten[{αβsolvalue[x], x}], {x, 0.4, x0}, PlotStyle -> {Thick, Red}], 
 AxesLabel -> {"α", "β", "x"}
]

enter image description here

Solution validity and error

To check that these solutions actually solve the equations I'll take some sample points

xvals = Range[0.001, 0.499, 0.001];
αβxvals = Flatten[{αβsol[#], {x -> #}}] & /@ xvals;

And @@ (eq1 /. αβxvals)
And @@ (eq2 /. αβxvals)

(* True
   True *)

Notice that I've left x == 0.5 out of these solutions, even though it technically lies within the range of the nontrivial solution. This is because as x -> x0 the solution branch converges to {α, β} == {0, 0}, which is undefined and causes the numerics to go a bit haywire (as shown in the next figure).

To get a better view of this, rearrange the equations once again to define:

f1[α_, β_, x_] := 470.44*(Coth[α] - 1/α) + 1329.74 (Coth[β] - 1/β) - α*x*840.92
f2[α_, β_, x_] := 664.87*(Coth[α] - 1/α) + 143.61 (Coth[β] - 1/β) - β*x*840.92

so that {α, β, x} solves the equations iff f1[α, β, x] == f1[α, β, x] == 0. Plotting these functions for {α, β} == αβsolvalue[x] shows that the error is indeed increasing as x -> x0:

GraphicsRow[{
  Plot[f1 @@ Flatten[{αβsolvalue[x], x}], {x, 0, x0}, 
   PlotRange -> All, PlotLabel -> "Error in equation 1"],
  Plot[f2 @@ Flatten[{αβsolvalue[x], x}], {x, 0, x0}, 
   PlotRange -> All, PlotLabel -> "Error in equation 2"]
  }, ImageSize -> 500]

enter image description here

So as x -> x0 the error in the solutions "explodes" to over 10^-11, which is large enough for Mathematica to decide that the values do not solve the equations.

The bifurcation at x == x0

Since we are looking for the zeros of f = {f1, f2} we can, with a bit of hand-waving, get some insight into the structure of the solution set by looking at the linearization of f = {f1, f2}.

f[α_, β_, x_] := {f1[α, β, x], f2[α, β, x]}
lin = D[f[α, β, x], {{α, β}}];

Evaluating in the limit as {α, β} -> {0, 0} and plotting the eigenvalues:

lin0 = Limit[lin, {α -> 0, β -> 0}]
Plot[Evaluate@Eigenvalues[lin0], {x, 0, 1}]

enter image description here

Clearly, one of the eigenvalues changes sign somewhere around x = 0.5, which indicates a change in the structure of the solutions (i.e.: a bifurcation). To pin down the bifurcation point, just solve for x:

x0 = x /. First@Solve[Eigenvalues[lin0][[2]] == 0, x]

(* 0.500003 *)

This is a long, long way from being a serious bifurcation analysis, so please don't quote me on it. Without more thought it grinds to a halt here because f is not actually defined at {α, β} == {0, 0}, so a whole bunch of tricks don't work and the analysis is on rather shaky theoretical grounds. However, something like this could work. And based on the numerics, it looks very much like there's a subcritical fold bifurcation at x0, and x0 does appear to give a pretty accurate upper x limit for the existence of solutions:

Plot[αβsolvalue[x], {x, x0 - 10^-6, x0}, AxesLabel -> {"x", "α,β"}]

enter image description here

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ With x=.5, FindRoot for alpha and beta around 0.1 will yield Alpha = 0.0113836 and Beta = 0.00677112. For Alpha and Beta around 0.00001, Alpha= -3.97303*10^-6 and Beta=-2.36328*10^-6. The plot clearly shows the intersection at alpha=beta=0 and the fact that the FindRoot result changes when you change the value you want the root around it, tells me that it is just an approximation done by the software. mathematically, however, I believe the only real result is alpha=beta=0. $\endgroup$ – Navid Rajil Jan 11 '18 at 6:21
  • 1
    $\begingroup$ @NavidRajil I believe the first plot in my answer suggests otherwise -- those two surfaces intersect, whether we have access to a closed form solution or just a numerical one. At those intersection points α and β are such that x1[α, β] == x2[α, β] == x, and so {x, α, β} must solve the original equations. $\endgroup$ – aardvark2012 Jan 11 '18 at 7:17
  • $\begingroup$ The equations you construct with x1[α, β] and x2[α, β] are actually mathematically undetermined $\endgroup$ – Navid Rajil Jan 11 '18 at 8:29
  • $\begingroup$ Aardvark2012, The equations you construct with x1[α, β] and x2[α, β] are actually mathematically undetermined at alpha and beta equal to 0. The reason is that (Coth[α]-1/α)/α and (Coth[β]-1/β)/α are 0/0 when α->0, so as the counter parts with α<->β. So any result that the machine gives, is just an approximation. the real mathematical answer is zero for both alpha and beta. You can notice it if you increase your precision and make it go closer to zero. $\endgroup$ – Navid Rajil Jan 11 '18 at 8:41
  • 2
    $\begingroup$ @NavidRajil Leaving (0,0) aside, if you check aardvark2012's nonzero solutions, you'll see they satisfy the OP's form of the equations, which shows there are at least two solutions. Given that (0,0) a removable discontinuity, hence special, my guess is that the nonzero solutions are more important. $\endgroup$ – Michael E2 Jan 11 '18 at 13:51
2
$\begingroup$

You can use FindRoot if you where know the approximate result would be,

FindRoot[{α*x*840.92 == 470.44*(Coth[α] - 1/α) + 
1329.74*(Coth[β] - 1/β), β*x*840.92 == 664.87*(Coth[α] - 1/α) + 143.61*(Coth[β] - 1/β)}, {{α, 1}, {β, 1}}]

the plot will tell you:

Plot3D[{\[Alpha]*x*840.92, 470.44*(Coth[\[Alpha]] - 1/\[Alpha]) + 1329.74*(Coth[\[Beta]] - 1/\[Beta]), 664.87*(Coth[\[Alpha]] - 1/\[Alpha]) + 143.61*(Coth[\[Beta]] - 1/\[Beta]), \[Beta]*x*840.92}, {\[Alpha], 0, 21}, {\[Beta], 0, 11}, PlotLegends -> {"\[Alpha]*x*840.92", "470.44*(Coth[\[Alpha]]-1/\[Alpha])+1329.74*(Coth[\[Beta]]-1/\\[Beta])", "664.87*(Coth[\[Alpha]]-1/\[Alpha])+143.61*(Coth[\[Beta]]-1/\[Beta]\)", "\[Beta]*x*840.92"}, AxesLabel -> Automatic]

My guess is both around zero or actually zero.

enter image description here

Manipulate[
Plot[{\[Alpha]*x*840.9, 
470.44*(Coth[\[Alpha]] - 1/\[Alpha]) + 
1329.74*(Coth[\[Beta]] - 1/\[Beta]), \[Beta]*x*840.92, 
664.87*(Coth[\[Alpha]] - 1/\[Alpha]) + 
143.61*(Coth[\[Beta]] - 1/\[Beta])}, {\[Alpha], -100, 
100}], {\[Beta], 0.0000001, 100}]

This Manipulate can let you change beta while you are plotting all RHS and LHS of your equations with respect to alpha. There is no place that they coincide except alpha=0 and beta=0. the other solution happens at {\[Alpha] -> 19.3646, \[Beta] -> 9.01653}. if you notice, orange and blue surfaces must intersect at those values and green and red also must intersect at those values. so, edited my answer and made it wider to have the other correct result mentioned.

| improve this answer | | | | |
$\endgroup$
  • $\begingroup$ Im assuming I dont know where the roots are. $\endgroup$ – Giovanni Baez Jan 10 '18 at 23:53
  • $\begingroup$ I was adding the approximation part using Plot3D. $\endgroup$ – Navid Rajil Jan 11 '18 at 0:54
  • $\begingroup$ I double checked, the final answer is alpha=0 and beta=0. $\endgroup$ – Navid Rajil Jan 11 '18 at 1:08
  • $\begingroup$ For these equations, alpha = 0 and beta = 0 are true and the values are independent of the value of x. This can be seen by Limit[Coth[x] - 1/x, x -> 0] = 0 making the rhs and lhs of each equation 0 when alpha,beta = 0. When using NSolve you generally do better with exact values in your equations, but in this case, FindRoot is the better choice. $\endgroup$ – Bill Watts Jan 11 '18 at 1:22
  • $\begingroup$ The Manipulate has a list of four components, yet shows only two curves. It might be instructive to understand the implications of that. $\endgroup$ – Daniel Lichtblau Jan 13 '18 at 15:46
0
$\begingroup$
eqns = α*x*840.92 == 
      470.44*(Coth[α] - 1/α) + 
       1329.74*(Coth[β] - 1/β) && β*x*840.92 == 
      664.87*(Coth[α] - 1/α) + 
       143.61*(Coth[β] - 1/β) // Rationalize // Simplify;

x = 1/10;

{#, FindRoot[eqns, {{α, 10^-9, 10^-6}, {β, 10^-9, 10^-6}}, 
     WorkingPrecision -> #] // N} & /@ Range[25, 125, 50]

(* {{25, {α -> -9.34083*10^-23, β -> -7.50174*10^-23}}, 
    {75, {α -> 3.73705*10^-55, β -> 3.00204*10^-55}}, 
    {125, {α -> -1.24106*10^-79, β -> 1.45138*10^-80}}} *)

The higher the precision, the closer the values are to zero. This indicates that the roots are zero.

| improve this answer | | | | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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