# Solving a transcendental equation with two parameters

I'm trying to solve an equation for T where T and two parameters M and S are all positive.

I tried solving it by evaluating the following expression, but it just gives me the same expression back. Any help would be much appreciated.

Solve [M*S + (T/0.755) - (T/0.755)*cosh[M/(2T/0.755)] == 0, T]

• First of all, cosh should be Cosh. Secondly, this is a transcendental equation that likely can't be solved analytically. You will probably need to use FindRoot with numbers in place of M and S. – march Jan 4 '18 at 4:51
• can you give an example on how to do that? @march – Jaffer Alfahdawi Jan 4 '18 at 4:53
• Perhaps the following would be useful. Some algebra produces: S = - (T/0.755)(1 - Cosh[2 T/0.755)/M, allowing S to be determined for any pair of {M,T}. Then Plot3D could be used to show the values {M,T,S} that satisfy the equation. e.g. Plot3D[-(T/0.755) ( 1 - Cosh[M/(2 T/0.755)] )/ M, {M, 0, 2}, {T, 0, 2}] – user6546 Jan 4 '18 at 5:15

One approach is to turn your equation into an ODE and then use NDSolve to to integrate.

eqn = M S+(T/0.755)-(T/0.755)*Cosh[M/(2T/0.755)]==0;


This can be simplified by replacing T with t M:

seqn = Assuming[M>0, eqn /. T->t M //Simplify];
seqn //TeXForm


$1. S+1.3245033112582782 t=1.3245033112582782 t \cosh \left(\frac{0.37749999999999995}{t}\right)$

This transcendental equation doesn't appear to be one that Mathematica can solve. However, we can turn the equation into an ODE and solve that:

ode = D[seqn /. t->t[S], S];
ode //TeXForm


$1.3245033112582782 t'(S)+1.=1.3245033112582782 t'(S) \cosh \left(\frac{0.37749999999999995}{t(S)}\right)-\frac{0.49999999999999994 t'(S) \sinh \left(\frac{0.37749999999999995}{t(S)}\right)}{t(S)}$

To use NDSolveValue we will also need an initial condition. For this we can use Solve:

t1 = t /. First @ Solve[seqn && 0<t<10 /. S->1, t]


Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.

0.153051

Now, we are ready to use NDSolveValue:

sol = NDSolveValue[{ode, t[1] == t1}, t, {S, 0, 10}];


NDSolveValue::ndsz: At S == 9.521862063632699*^-7, step size is effectively zero; singularity or stiff system suspected.

Here is a plot:

Plot[sol[t], {t, 0, 10}]


The corresponding solution for T is:

Tsol[s_, m_] := m sol[s]


Here is a Plot3D of the LHS of the equation showing that it is satisfied to less than 10^-4. Using exact inputs and higher accuracy/precision goals can produce a better result:

Plot3D[eqn[[1]] /. T->Tsol[S, M], {S, 0, 10}, {M, 0, 10}]


• So what would be the formula for T if I want different values for M and S? @Carl Woll – Jaffer Alfahdawi Jan 4 '18 at 5:45
• The reason I want formula for T is because I want to use it in a android application to calculate value for T given values for M and S. Cheers – Jaffer Alfahdawi Jan 4 '18 at 5:48
• @JafferAlfahdawi You can look at questions like this for help in using NDSolve output outside of Mathematica. – Carl Woll Jan 4 '18 at 5:54

You can get a symbolic formula approximating the solution via Chebyshev series, which can be computed as shown in About multi-root search in Mathematica for transcendental equations.

As @Carl Woll observed, your equation is a function of m == M/T and S. If we solve for m as a function of S, we can then get a formula for T = M/m.

From a comment, we get the domain for S needed to compute the Chebyshev series:

M can be between 1 and 200 and S can be between 0.01 to 0.07.

Below, we used FourierDCT to compute the Chebyshev coefficients. We have to find the values of m as S ranges over the Chebyshev nodes. We use FindRoot to do this with high precision, so that we can construct a near machine-precision approximation to m, which is given by m = mFN[S].

eqn = M S + (T/0.755) - (T/0.755)*Cosh[M/(2 T/0.755)] // Rationalize;
eqn3 = eqn/T /. {M -> m T} // Simplify; (* m == M/T, or T == M/m *)

Block[{n = 32, S1, S2, S, snodes, (*coeffs,*) coeffs2, mvals, m},
{S1, S2} = Rationalize@{0.01, 0.07};  (* domain for S *)
snodes = Rescale[Sin[Pi/2 Range[n, -n, -2]/n], {-1, 1}, {S1, S2}];
mvals = m /. Rest@FoldList[
FindRoot[eqn3 /. S -> #2, {m, m /. #1, 1/10, 1000}, WorkingPrecision -> 40] &,
{m -> 10 Sqrt[S1]}, snodes];
coeffs2 = Sqrt[2/n] FourierDCT[mvals, 1];
coeffs2[[{1, -1}]] /= 2;
coeffs = Drop[coeffs2,  (* drop coefficients smaller than machine precision *)
Module[{sum = 0}, -LengthWhile[
Reverse@coeffs2, (sum += Abs[#]) < Abs[First@coeffs2] $MachineEpsilon &]]]; With[{s1 = S1, s2 = S2, cc = coeffs, k = Length@coeffs - 1}, mFN = Function[{S}, (* solution for m *) cc.Cos[Range[0, k] ArcCos[Rescale[S, {s1, s2}, {-1, 1}]]] ]] ]; tFN = Function[{M, S}, M/mFN[S]]; (* solution for t *)  Answer: Here is the double-precision formula for t that is an approximate solution to the OP's equation: N[tFN[M, S]] (* formula for t *) (* M/(0.4221956380745171 + 0.31559498734947056 (-1.3333333333333333 + 33.333333333333336 S) - 0.0007452700193212066 Cos[2. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] - 0.00008989046047576488 Cos[3. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] + 9.349727403445521*^-7 Cos[4. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] + 6.316071780569096*^-8 Cos[5. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] - 1.2738992410763405*^-9 Cos[6. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] - 5.477086714861503*^-11 Cos[7. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] + 1.7873875250882988*^-12 Cos[8. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] + 4.9699383243064014*^-14 Cos[9. ArcCos[-1.3333333333333333 + 33.333333333333336 S]] - 2.529762133661728*^-15 Cos[10. ArcCos[-1.3333333333333333 + 33.333333333333336 S]])  Simple check of accuracy: is there a small residual? eqn /. T -> tFN[M, S] /. {M -> RandomReal[{1, 200}, 1000], S -> RandomReal[{0.01, 0.07}, 1000]} // Abs // Max (* 4.54747*10^-13 *)  Somewhat more sophisticated look: We compare the residual error with the largest of the three terms in eqn, divided by $MachineEpsilon:

Total[#, {2}]/Max /@ Abs[#] &@Transpose[
List @@ eqn /.
T -> tFN[M, S] /. {M -> RandomReal[{1, 200}, 10000],
S -> RandomReal[{0.01, 0.07}, 10000]}
]/\$MachineEpsilon //
Histogram[#, Automatic, "LogCount",
AxesLabel -> {"Error (ϵ)", "Count"},
PlotLabel -> "Error distibution relative to largest term of eqn ε_mach"] &


Almost half of the values, 4600 out of 10000, have zero error. Most of the rest have a 1 ulp error.

Update: avoiding transcendental functions

I forgot to mention Clenshaw's recurrence for evaluating a Chebyshev series. The code used in mFN is an efficient way to compute mFN[S] in the Mathematica computational environment, even with the trigonometric functions. However, it is just a polynomial and can be computed without a math library. The best way is to use the Clenshaw recurrence (see chebeval here for instance).

If desired, one can see the polynomial expression by changing the code for mFM to the following:

mFN = Function[{S}, cc.ChebyshevT[Range[0, k], Rescale[S, {s1, s2}, {-1, 1}]]]


If you divide your equation by T you get an equation which can easily solved for S:

gl = M S + (T/0.755) - (T/0.755)*Cosh[M/(2 T/0.755)];
glT = Simplify[gl /T /. M -> mdT T ]
s[mdT_] := Evaluate[S /. Solve[glT == 0, S][[1]]]
Plot[s[mdT ], {mdT, -10, 10}, AxesLabel -> {m/T, S}](-1.3245 + 1.3245 Cosh[0.3775 M/T])/(M/T)


S is a function of M/T. There is no need to solve for T.

If you know the range of M,S you can use InverseSeries[] to get an approximation of the InverseFunction:

M/Normal[ InverseSeries[Series[s[mdT], {mdT, 0, 3}], S] ](* ~T *)
(* T=M/(10.596 S - 14.128 S^3) *)
`
• M can be between 1 and 200 and S can be between 0.01 to 0.07. M is the span length of conductor between two poles and S is the conductor sag percentage and T is the tension of the conductor if that helps. – Jaffer Alfahdawi Jan 4 '18 at 14:23