1
$\begingroup$

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]
Improve this question
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – march
    Commented Jan 4, 2018 at 4:51
  • $\begingroup$ can you give an example on how to do that? @march $\endgroup$
    – Jaffer Alfahdawi
    Commented Jan 4, 2018 at 4:53
  • $\begingroup$ 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}] $\endgroup$
    – user6546
    Commented Jan 4, 2018 at 5:15

3 Answers 3

Reset to default
4
$\begingroup$

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

Here is your equation:

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

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ So what would be the formula for T if I want different values for M and S? @Carl Woll $\endgroup$
    – Jaffer Alfahdawi
    Commented Jan 4, 2018 at 5:45
  • $\begingroup$ 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 $\endgroup$
    – Jaffer Alfahdawi
    Commented Jan 4, 2018 at 5:48
  • $\begingroup$ @JafferAlfahdawi You can look at questions like this for help in using NDSolve output outside of Mathematica. $\endgroup$
    – Carl Woll
    Commented Jan 4, 2018 at 5:54
3
$\begingroup$

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

Mathematica graphics

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}]]]
Improve this answer
$\endgroup$
1
$\begingroup$

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.

enter image description here

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) *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Jaffer Alfahdawi
    Commented Jan 4, 2018 at 14:23

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.