# Solve a transcendental equation in Mathematica

I believe it is called a "transcendental equation", and I would like to find the best possible solution.

I tried to use Solve:

Solve[Cos[ϕ0]/ϕ0 == h/(2 s0), ϕ0]


During evaluation of In[5]:= Solve::nsmet: This system cannot be solved with the methods available to Solve. >>

Solve[Cos[ϕ0]/ϕ0 == h/(2 s0), ϕ0]


And it didn't work. Since this is the first time I am working with this kind of equation, could somebody help me? I know I can use Taylor expansion for h/s0 << 1 or h/s0 >> 1, but I don't want to do that. The reason why not is because this would not be a realistic model. In my case h is approximately the same as s0.

### Edit

BTW, I would be satisfied with numerical solutions as well -- but only in case that is the only option.

• Numerics is your only option, yes… – J. M. will be back soon Jul 8 '15 at 8:33
• Ok. If so.. How? I though NSolve would do the job, but it doesn't. Problem: "This system cannot be solved with the methods available to NSolve. " – skrat Jul 8 '15 at 8:43
• Use FindRoot[], and have a good initial guess ready. – J. M. will be back soon Jul 8 '15 at 8:44
• For any numerical value of the right-hand side, Solve does a nice job when you add the domain Reals, such as Solve[Cos[ϕ0]/ϕ0 == 1/4, ϕ0, Reals] – Fred Simons Jul 8 '15 at 10:13
• This message is not an error message. It is information for the user how the solution is obtained: first the inexact (machine) numbers are replaced with exact numbers, then Solve finds the exact solution of the resulting equation and then returns the numerical approximation of the exact result. Keep in mind that the solution set of an equation can be highly senitive for small changes in the coefficients. – Fred Simons Jul 15 '15 at 11:26

You're in luck if you're only interested in a domain where x = h/(2 s0) is approximately 1/2. That domain is contained by the domain of InverseFunction applied to Cos[φ0]/φ0:

ϕ = InverseFunction[Cos[#]/# &]


InverseFunction can only work on a domain on which the function is one-to-one. In this case, it's from zero to the first minimum around ϕ0 -> 2.79839:

Plot[Cos[ϕ0]/ϕ0, {ϕ0, 0, 10}]


Now the equation has a unique solution for h/(2 s0) greater than the first maximum, which about h/(2 s0) == 0.161228 or for 0 < ϕ0 < 1.35119:

{minhs, foo} = FindMaximum[Cos[ϕ0]/ϕ0, {ϕ0, 6}]
{maxϕ} = NSolve[Cos[ϕ0]/ϕ0 == minhs && 0 < ϕ0 < 2, {ϕ0}]
(*
{0.161228, {ϕ0 -> 6.12125}}
{{ϕ0 -> 1.35119}}
*)


If these limitations are acceptable, then the inverse function is your solution:

ϕ0 -> ϕ[h/(2 s0)]


It is a symbolic solution and computable to arbitrary precision.

ϕ[0.5]
(*  1.02987  *)

N[ϕ[1/2], 20]
(*  1.0298665293222588276  *)

Plot[ϕ[x], {x, 0, 1},
GridLines -> {{minhs}, {ϕ0 /. maxϕ}}, Mesh -> {{minhs}},


The lower right region shows the valid solution.

Alexei Boulbitch's answer is quite nice, but here are another couple of useful techniques for situations like this. First, if all you need is a graph of the function, you can use ParametricPlot to get that. You want $\phi(x)$, where $\cos \phi/\phi = x$; so you can use $\phi$ as a parameter and plot the curve $(\cos \phi/\phi, \phi)$:

ParametricPlot[{Cos[ph]/ph, ph}, {ph, 0, Pi/2}, PlotRange -> {{0, 1}, Automatic}, AspectRatio -> 1/GoldenRatio]


Second, you can use Interpolation to return an approximate inverse function. Generate a list of points you want your graph to go through; this will be a table of points of the form $(\cos \phi/\phi, \phi)$, as above. Then feed it into Interpolation:

pts = Table[{Cos[ph]/ph, ph}, {ph, Pi/40, Pi/2, Pi/40}];
phiroot[x_] = Interpolation[pts][x];


phiroot now returns the function $\phi(x)$ that you were looking for:

N[phiroot[1/2]]

(* 1.02986 *)


For comparison, the value found by FindRoot differs from this value by about one part in $10^5$:

FindRoot[Cos[x]/x == 1/2, {x, 1}]

(* {x -> 1.02987} *)


The weakness of this approach is that if you try to find phiroot for a value of x outside the range of pts, then you'll be extrapolating rather than interpolating and your results will go pear-shaped rather quickly. Still, if you're only concerned with a fixed range of $x$, this approach works rather well.

Looking at the plot:

Manipulate[Plot[Cos[x]/x - y, {x, 0, 2 \[Pi]}], {y, 0, 1}]


that is, here:

and varying the parameter y=h0/(2s0), one finds that the first root is close to x==1. Let us look for this root only, than it is easy to find a list of values with the structure {y,x0}, where y is the value of the parameter and x0 is the root:

    lst = Table[{y, FindRoot[Cos[x]/x - y == 0, {x, 1}][[1, 2]]}, {y, 0,
1, 0.05}]

(*  {{0., 1.5708}, {0.05, 1.49593}, {0.1, 1.42755}, {0.15, 1.36464}, {0.2,
1.30644}, {0.25, 1.25235}, {0.3, 1.20191}, {0.35, 1.15474}, {0.4,
1.11051}, {0.45, 1.06896}, {0.5, 1.02987}, {0.55, 0.993021}, {0.6,
0.958252}, {0.65, 0.925404}, {0.7, 0.894337}, {0.75,
0.864927}, {0.8, 0.837061}, {0.85, 0.810634}, {0.9,
0.785553}, {0.95, 0.76173}, {1., 0.739085}}  *)


Then this:

    model = a/(1 + b*y);
ff = FindFit[lst, model, {a, b}, y]
Show[{
ListPlot[lst],
Plot[model /. ff, {y, 0, 1}, PlotStyle -> Red]
}]


returns the fitting parameters for the model:

(* {a -> 1.58663, b -> 1.10311}  *)


and a plot for the visual inspection of the fitting quality:

Though it is approximate, the obtained solution:

x=1.59/(1+1.1*y)

is analytic, simple and rather accurate. The other solution can be found similarly.

Have fun!

• I have to ask you one question though... Where did you get that model? Based on what? Why didn't you for example go for a third degree polynomial or something..? – skrat Jul 8 '15 at 11:17
• Also, please help me understand what this does. I will walk you through and let me know if my understanding is wrong. Basically what you did is you firstly created a list based on numerical solution of the original equation for different parameter values y. Than you fitted a function which fits the best and the end result is a nice analytical solution (which is not exact but good) to the original equation? Right? – skrat Jul 8 '15 at 11:31
• @skrat Right. "Where did you get that model?" Just proposed one that I expected to behave like what I saw in the output of the ListPlot. You could also propose another one. An infinite number of models will fit this list. In this situation, and since the model is anyway an approximation, the simplest one is the best. – Alexei Boulbitch Jul 8 '15 at 12:29