Finding the analytic root of a multivariate function that involves trigonometric functions

Question

I am currently stuck in a problem, where I need to find an analytic expression $x_0=f(\alpha,\beta)$. Let me first try to explain what $x_0$ is:

Let $x$ be a real-valued variable, which must obey the following relation: \begin{equation} x^2 + \left[a_1(\omega) \alpha + a_2(\omega) \beta \right] x + \left[ a_3(\omega) + a_4(\omega) \alpha \beta \right] \leq 0 \end{equation} with $\omega \in [0,\pi]$ and both $\alpha$ and $\beta$ are positive and real-valued. The $a_i(\omega)$ terms are defined as follows: \begin{equation} a_1(\omega) = \frac{4 \left( cos(\omega) - 1\right)}{sin^2(\omega) + 2 \left(\cos(\omega) - 1\right)} \end{equation} \begin{equation} a_2(\omega) = \frac{-8}{cos(2\omega)-1} \end{equation} \begin{equation} a_3(\omega) = \frac{3-cos(\omega)}{cos(\omega)-1} \end{equation} \begin{equation} a_4(\omega) = \frac{-16}{ sin^2(\omega) \left(cos(\omega)-1\right)} \end{equation} The equation for $x$ can be written in terms of \begin{equation} x^2 + p x + q \leq 0 \end{equation} with \begin{equation} p=a_1(\omega)\,\alpha + a_2(\omega)\,\beta \end{equation} \begin{equation} q=a_3(\omega) + a_4(\omega)\,\alpha\,\beta. \end{equation}

If $4 q < p^2$ is true, one can obtain a possible solution for $x$: \begin{equation} \frac{1}{2}\left(-\sqrt{p^2-4q} - p\right) \leq x \leq \frac{1}{2}\left(\sqrt{p^2-4q} - p\right) \end{equation}

Now, I am especially interested in finding the maximum of the upper bound of this solution regarding $\omega$, which I call $x_0(\alpha,\beta)$: \begin{equation} x_0(\alpha,\beta) = \max_{\omega} \frac{1}{2}\left(\sqrt{p^2-4q} - p\right) \end{equation}

To estimate an expression for $x_0(\alpha,\beta)$, I thought I could calculate the first derivative of $(\sqrt{p^2-4q}-p)/2$, setting it zero and solve for $\omega$. This should give me an expression like $\omega_0 = f(\alpha,\beta)$, where $\omega_0$ is the $\omega$-value at which the maximum $x_0$ is located. This expression could then be inserted into $(\sqrt{p^2-4q}-p)/2$ and should result in an expression for $x_0=f(\alpha,\beta)$.

Of course, it should also be checked, if $x_0$ is actually a maximum via the second derivative. However, I checked the shape of $(\sqrt{p^2-4q}-p)/2$ for some $\alpha$ and $\beta$ values I am interested in and it didn't seem to have local minimas in that range.

Basically, my problem is the step, where I want to solve for $\omega$ when setting the first derivative to zero: \begin{equation} \frac{\partial}{\partial \omega} \frac{1}{2}\left(\sqrt{p^2-4q} - p\right) = 0 \end{equation} My problem is, that $\omega$ occurs only in the trigonometric functions in the $a_i(\omega)$-terms. I tried to solve the problem with with a symbolic programming library (sympy), but failed at this point to estimate a solution.

Has anyone a hint for me, how such a problem can be solved or in which direction I could go? Maybe there is some kind of trick I can use to simplify the problem?

Thank you all in advance!

Answer 1

I haven't double checked to make sure there are no typing errors but something like that will do the trick. The main idea was to write a polynomial system in order to compute the Groebner basis in lexicographic order.

from sympy import *

a, b, w = symbols("a b w")

a1 = (4*(cos(w)-1))/(sin(w)**2+2*(cos(w)-1))

a2 = (4)/(sin(w)**2)

a3 = (3-cos(w))/(cos(w)-1)

a4 = (-16)/(sin(w)**2*(cos(w)-1))

p = a1 * a + a2 * b

q = a3 + a4 * a * b

x, r, s = symbols("x r s")

EQN_1 = cancel(x**2 + x*p + q)

EQN_2 = cancel(diff(q,w) + x*diff(p,w))

num_EQN_1, den_EQN_1 = fraction(EQN_1)

num_EQN_2, den_EQN_2 = fraction(EQN_2)

EQNS = Matrix([[num_EQN_1, num_EQN_2, r**2+s**2-1]])
EQNS = EQNS.subs(sin(w),r)
EQNS = EQNS.subs(cos(w),s)

"""
SOL = groebner(EQNS, s, r, x, order='lex')
SOL = groebner(EQNS, x, s, r, order='lex')
"""
SOL = groebner(EQNS, x, r, s, order='lex')

"""
TESTING for alpha = beta = 1
"""

SOL_1 = SOL
SOL_1 = SOL_1.subs(a,1)
SOL_1 = SOL_1.subs(b,1)
sym_cos_w = solve(SOL_1[2],s)
num_cos_w = sym_cos_w

for i in range(0,len(sym_cos_w)):
    num_cos_w[i] = N(num_cos_w[i])
    if ( abs(im(num_cos_w[i])) < 1e-6 ) and ( abs(re(num_cos_w[i])) < ( 1 + 1e-6 ) ):
        num_sin_w = solve(SOL_1[1].subs(s,num_cos_w[i]),r)
        for j in range(0,len(num_sin_w)):
            num_x = SOL_1[0].subs(s,num_cos_w[i])
            num_x = solve(num_x.subs(r,num_sin_w[j]))
            for k in range(0,len(num_x)):
                print('cos_w = ', num_cos_w[i], 'sin_w = ', num_sin_w[j], 'x = ', num_x[k])

EDIT:

The answer was updated because the line EQN_2 = cancel(2*diff(q,w) + x*diff(p,w)) was wrong. The correct assignment is EQN_2 = cancel(diff(q,w) + x*diff(p,w)).

Thank you, FloSe, for pointing it out.