5
$\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!

Improve this question
$\endgroup$
4
  • 2
    $\begingroup$ Try FindRoot[{A1, B2, C3}, {{h, 1}, {L, 1}, {Umax, 1}}. $\endgroup$ Commented Mar 7, 2016 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$ Commented Mar 7, 2016 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$ Commented Mar 7, 2016 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$ Commented Mar 7, 2016 at 3:15

2 Answers 2

Reset to default
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.

Improve this answer
$\endgroup$
5
  • $\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$ Commented Mar 7, 2016 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$ Commented Mar 7, 2016 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$ Commented Mar 7, 2016 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$ Commented Mar 7, 2016 at 12:46
  • $\begingroup$ @user3185748 You are welcome, I'm glad you find my answer helpful. $\endgroup$ Commented Mar 7, 2016 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

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.