# Solution for coupled equations with trigonometric functions

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.

• 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. Jan 11 '18 at 2:14

## 3 Answers

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


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


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


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


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


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", "α,β"}]


• 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. Jan 11 '18 at 6:21
• @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. Jan 11 '18 at 7:17
• The equations you construct with x1[α, β] and x2[α, β] are actually mathematically undetermined Jan 11 '18 at 8:29
• 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. Jan 11 '18 at 8:41
• @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. Jan 11 '18 at 13:51

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.

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.

• Im assuming I dont know where the roots are. Jan 10 '18 at 23:53
• I was adding the approximation part using Plot3D. Jan 11 '18 at 0:54
• I double checked, the final answer is alpha=0 and beta=0. Jan 11 '18 at 1:08
• 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. Jan 11 '18 at 1:22
• The Manipulate has a list of four components, yet shows only two curves. It might be instructive to understand the implications of that. Jan 13 '18 at 15:46
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.