How to efficiently optimize a function with two local minima?

I want to minimize a function that I know it has two local minima, but using scipy to optimize it I don't get the global minimum.
I want to know of a way to do it that is fast, as I have to do it thousands of times.

I have functions to compute a result from two vectors of 3 components (v1 and v2) and an angle theta:

def get_v_theta(v, theta):
    vx = v[0]
    vy = v[1]
    vz = v[2]
    v_theta = (vx + vy)/2 + (vx - vy)/2*np.cos(2*theta) + vz*np.sin(2*theta)
    return v_theta


def get_result(v1_theta, v2_theta):
    v_max = max(v1_theta, v2_theta)
    v_min = min(v1_theta, v2_theta)
    v_med = (v_max + v_min)/2
    v_alt = (v_max - v_min)/2
    if v_med > 0:
        result = (v_alt/100 + v_med/400)**(-1)
    else:
        result = 100/v_alt
    return result


def get_result_theta(theta, v1, v2):
    v1_theta = get_v_theta(v1, theta)
    v2_theta = get_v_theta(v2, theta)
    result = get_result(v1_theta, v2_theta)
    return result

I want to find the angle theta that minimizes result, given two vectors.
For example, if my vectors are:

import numpy as np

v1 = np.array([2.91931, 0.758059, 0.10201])
v2 = np.array([-1.29405, 7.76101, -1.75803])

I can plot result vs theta and see that the theta that minimizes result should be around -1.4.

import matplotlib.pyplot as plt

thetas = np.arange(-np.pi/2, np.pi/2, np.pi/100)
result = []
for theta in thetas:
    result.append(get_result_theta(theta, v1, v2))
plt.plot(thetas, result)

But if I use scipy to find that value, it returns 0.147.

from scipy import optimize

res = optimize.minimize_scalar(get_result_theta, args=(v1, v2), method='bounded', bounds=(-np.pi/2, np.pi/2))

How can I get the global minimum instead of a local one, without taking much longer (I have to do this computation thousands of times), knowing that there should be only 2 local minima?

You can solve this problem analytically, and there are some nice simplifications of the equations. First, for ease of the discussion, instead of minimizing the result, you can maximize the inverse.

Now note that for any p, q we have max(p, q) + min(p, q) = p + q and max(p, q) - min(p, q) = |p - q|.

The functions V(Θ) are of the form linear-trigonometric

A cos(2Θ) + B sin(2Θ) + C

So Vmed is of that form and Valt is piecewise linear-trigonometric, with the the coefficients changing sign where

(A' - A") cos(2Θ) + (B' - B") sin(2Θ) + (C' - C") = 0.

There are explicit formulas to obtain these Θ, as well as the zeroes of Vmed (giving max(Vmed, 0), another piecewise LT function).

Finally, by combining the intervals you get the inverse of result as a piecewise linear-trigonometric function, and you find the stationary points by canceling the first derivative in every piece and keeping the maximum maximorum.