Solve a transcendental equation in Mathematica

Question

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.

Answer 1

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}}, 
 MeshShading -> {Dashed, Automatic}]

The lower right region shows the valid solution.

Answer 2

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.

Answer 3

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!