4
$\begingroup$

I'm currently trying to use Mathematica to solve a problem for my fluid mechanics class. I need to solve the following system of 3 equations:

$$\begin{align*} \dfrac{h}{a} &= \cot \left(\dfrac{h/a}{2m\,/(U_{\infty}a)} \right)\\ \dfrac{L}{a} &= \sqrt{1+\dfrac{2m}{U_{\infty}a}}\\ \dfrac{U_{\mathrm{max}}}{U_{\infty}} &= 1+\dfrac{2m\,/(U_{\infty}a)}{1+\frac{h^2}{a^2}} \end{align*}$$

The following is known:

$$ \dfrac{m}{U_{\infty}\,a}=1 $$

I am trying to find the following:

$$\begin{align*} \dfrac{h}{a} \\\dfrac{L}{a} \\\dfrac{L}{h} \\\dfrac{U_{\mathrm{max}}}{U_{\infty}} \end{align*}$$

I have inserted the equations into Mathematica:

A1 = h/a == Cot[(h/a)/((2*m)/(U*a))]
B2 = L/a == Sqrt[1 + ((2*m)/(U*a))]
C3 = Umax/U == 1 + (((2*m)/(U*a))/(1 + ((h^2)/(a^2))))

I then set some initial conditions, these should not matter so long as their ratio is 1:

m = 1
a = 1
U = 1

I now have a system of 3 equations with 3 unknowns so I tried solving with 2 methods but I get:

Solve[{A1, B2, C3}, {h, L, Umax}, Reals]

Gives me:

This system cannot be solved with the methods available to Solve.

Also:

FindInstance[{A1, B2, C3}, {h, L, Umax}, Reals]

Gives me:

The methods available to FindInstance are insufficient to find the requested instances or prove they do not exist. 

I know the answers should give me:

$$\begin{align*} \dfrac{h}{a} &= 1.307\\ \dfrac{L}{a} &= 1.732\\ \dfrac{L}{h} &= 1.326\\ \dfrac{U_{\mathrm{max}}}{U_{\infty}} &= 1.739 \end{align*}$$

How else could I try solving this problem to give me the ratios I need?

Thanks!

$\endgroup$
  • 2
    $\begingroup$ Try FindRoot[{A1, B2, C3}, {{h, 1}, {L, 1}, {Umax, 1}}. $\endgroup$ – march Mar 7 '16 at 3:03
  • $\begingroup$ Sure enough that did it! Thanks so much! If you care to create a solution I'll accept it. (Note: missing closing ] for FindRoot[]) $\endgroup$ – user3185748 Mar 7 '16 at 3:08
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Mar 7 '16 at 3:09
  • $\begingroup$ In general, transcendental equations don't do so well in Solve[]; FindRoot[] with a good starting value does better in those cases. $\endgroup$ – J. M. is away Mar 7 '16 at 3:15
7
$\begingroup$

There is a smart use of Solve which eliminates certain variables and solves the system with respect to the remaining ones, see e.g. this answer : Efficient code for solve this equation.
However, here it appears to take too much time. Moreover such an approach to this system cannot be quite successful, since there are infinitely many roots of certain transcendetal equation present in the system.

Nevertheless we can eliminate most variables using Eliminate.

Eliminate[{ h/a == Cot[h/a/(2 m/(U a))], 
            L/a == Sqrt[1 + 2 m/(U a)], 
            Um/U == 1 + (2 m)/(U a)/(1 + h^2/a^2), 
            m/(U a) == 1, 
            A1 == h/a, 
            B2 == L/a, 
            C3 == Um/U},
           {a, h, m, U, Um}, InverseFunctions -> True]
   B2^2 == 3 && (1 + A1^2) C3 == 3 + A1^2 && 2 ArcCot[A1] == A1

Now we can demonstrate how to use an exact method like Solve. The full system reduces to another one much simpler:

{ A1 == Cot[A1/2],  (1 + A1^2) C3 == 3 + A1^2 }

Nevertheless we couldn't get the solutions this way

Solve[{A1 == Cot[A1/2],   (1 + A1^2) C3 == 3 + A1^2}, {A1, C3}, Reals]
Solve::nsmet: This system cannot be solved with the methods available to Solve. >>

since there are infinitely many solutions of this simple (however a transcendental one) equation A1 == Cot[A1/2]. If we restrict A1 to a given range e.g. 0 < A1 < 6 Pi, we can find all corresponding solutions (there are three of them represented in terms of the Root objects):

pts = {A1, Cot[A1/2]} /. Solve[{A1 == Cot[A1/2], 0 < A1 < 6 Pi}, A1]
Plot[{A1, Cot[A1/2]}, {A1, 0, 6 Pi}, Epilog -> {Red, PointSize[0.01], Point[pts]}]

enter image description here

These are exact solutions ( see e.g. How do I work with Root objects?).

Warning: I found a puzzling bug (in Mathematica 10.1 as well as in 10.3): compare e.g. A1 /. Solve[{A1 == Cot[A1/2], 0 < A1 < 799/100 Pi}, A1] or any other number close to it with A1 /. Solve[{A1 == Cot[A1/2], 0 < A1 < 8 Pi}, A1]. The latter finds only one root while it should find four roots.

From the context of the problem we are interested in finding the root in the range 0 < A1 < Pi/2, thus we supplement the system with such a condition:

{A1, C3} /. 
Solve[(1 + A1^2) C3 == 3 + A1^2 && 2 ArcCot[A1] == A1 && 0 < A1 < Pi/2,
      {A1, C3}]

N @ %
  {{ Root[{-2 ArcCot[#1] + #1 &, 1.30654237418880620223}], 
      (3 + Root[{-2 ArcCot[#1] + #1 &, 1.30654237418880620223}]^2)/(
        1 + Root[{-2 ArcCot[#1] + #1 &, 1.30654237418880620223}]^2)}}

     {{1.30654, 1.73881}}

and since B2 == Sqrt[3] i.e. Sqrt[3] // N yields 1.73205.

$\endgroup$
  • $\begingroup$ I now understand why you wanted a transcentals tag...I merely sought to simplify and search for root near sol. Poor I guess...your answer is instructive so +1 $\endgroup$ – ubpdqn Mar 7 '16 at 10:38
  • $\begingroup$ @ubpdqn Thanks. Recently I rarely answer the questions, but sometimes I still find certain reasons. Here it was the inability of Solve to eliminate most of variables even though I used Solve[{eqs}, {vars}, {varsToElim}] where I supplemented eqs with appropriate conditions. $\endgroup$ – Artes Mar 7 '16 at 10:55
  • $\begingroup$ I think your answer and link show how to use Solve and points to Root objects and is better than just getting the desired answer which was my lazy contribution. Understanding trumps lucky special cases $\endgroup$ – ubpdqn Mar 7 '16 at 11:00
  • $\begingroup$ Wow, thanks so much for the comprehensive answer. You offered a number of ideas that certainly hadn't occurred to me. Thanks so much! $\endgroup$ – user3185748 Mar 7 '16 at 12:46
  • $\begingroup$ @user3185748 You are welcome, I'm glad you find my answer helpful. $\endgroup$ – Artes Mar 7 '16 at 12:55
3
$\begingroup$

Using what is given and simplifying:

la = N@Sqrt[3]
ha = x /. FindRoot[Cot[x/2] - x, {x, 1}]
u = 1 + 2/(1 + ha^2)
lh = la/ha

yielding

Grid[{la, ha, u, lh} /. {la -> {"l/a", la}, ha -> {"h/a", ha}, 
   u -> {"Umax/U", u}, lh -> {"l/h", lh}}]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.