Wrong solution to a quadratic equation

Question

I'm facing a mathematical/programming problem which I don't really understand. In the following α,β,θc are given parameters, ϕ is the variable. I have to solve the equation Cp=0, where:

and:

The solutions that I have computed are:

Coding the above solutions with Python/Numpy, I find that:

If β=0, then both solutions are correct: Cp(ϕ1)=Cp(ϕ2)=0.
However, if β!=0 then Cp(ϕ1)!=0 (which is wrong) and Cp(ϕ2)=0. Why? I checked again and again the signs and terms on my coded expressions, they seem fine to me...

Here is the code:

import numpy as np
theta_c = np.deg2rad(10)
alpha = np.deg2rad([1, 2, 4, 6, 8, 10, 12, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85])
beta = np.deg2rad(15)

lamb = np.cos(alpha) * np.cos(beta)
nu = -np.sin(beta)
omega = np.sin(alpha) * np.cos(beta)

common_1 = lamb * nu * np.tan(theta_c)
common_2 = omega * np.sqrt(omega**2 + nu**2 - lamb**2 * np.tan(theta_c)**2)
common_3 = omega**2 + nu**2

phi_1 = np.arcsin((-common_1 + common_2) / common_3)
phi_2 = np.arcsin((-common_1 - common_2) / common_3)

Cp = lambda phi: lamb * np.sin(theta_c) + nu * np.cos(theta_c) * np.sin(phi) - omega * np.cos(phi) * np.cos(theta_c)

print(Cp(phi_1))
# OUTPUT: this is wrong!
# [-0.02357418 -0.04432021 -0.0779439  -0.10278244 -0.12111675 -0.13494261
#  -0.1457658  -0.15858164 -0.17562354 -0.19112922 -0.20703468 -0.22411243
#  -0.24265822 -0.26274848 -0.28434621 -0.3073487  -0.33161092 -0.35695785
#  -0.38319138 -0.41009444 -0.43743318 -0.4649581 ]

print(Cp(phi_2))
# OUTPUT: all values are very close to zero. Good solution.
# [-1.04083409e-17 -4.16333634e-17  0.00000000e+00 -1.38777878e-17
#  -2.77555756e-17 -2.77555756e-17  0.00000000e+00  0.00000000e+00
#   5.55111512e-17 -1.11022302e-16  1.11022302e-16  5.55111512e-17
#  -1.11022302e-16  0.00000000e+00 -3.88578059e-16 -1.11022302e-16
#  -2.77555756e-16  1.66533454e-16  5.55111512e-17  2.22044605e-16
#   5.55111512e-17 -2.22044605e-16]

Alternatively, I could use fsolve from Scipy optimize, which I have already verified it computes the correct solutions. However, this procedure will be executed many times by another iterative algorithm, so I'd rather use the analytical solutions.

user_17892216 · Answer

Seems like your analytical solution has a mistake.

I rewrote the target function in form

Cp(ϕ) = 2 * (A + B sin ϕ - C cos ϕ)^2

Solutions of Cp = 0 come from A + B sin ϕ - C cos ϕ = 0 equation, which I solved using Weierstrass substitution (WS). The solutions ϕ₁, ϕ₂ I found are

# t is parameter in WS

            -B ± sqrt(B^2 - (A^2 - C^2))
t₁ and t₂ = ----------------------------
                      A + C

ϕ₁ = 2atan(t₁)
ϕ₂ = 2atan(t₂)

Here is code snippet written in Julia

using Plots

θ = deg2rad(10)
α = deg2rad(10)
β = deg2rad(15)

λ = cos(α) * cos(β)
ν = -sin(β)
ω = sin(α) * cos(β)

# Cp = 2 * (A + B sin ϕ - C cos ϕ)^2
A = λ * sin(θ)
B = ν * cos(θ)
C = ω * cos(θ)

t₁ = (-B + sqrt(B^2 - (A^2 - C^2)))/(A + C)
t₂ = (-B - sqrt(B^2 - (A^2 - C^2)))/(A + C)

ϕ₁ = 2atan(t₁)
ϕ₂ = 2atan(t₂)

Cfunc(ϕ) = 2*(A + B * sin(ϕ) - C * cos(ϕ))^2

plot(Cfunc; xlims=(-10, 10), label="Cₚ")
scatter!([ϕ₁, ϕ₂], Cfunc.([ϕ₁, ϕ₂]); label="Solutions")

producing

I played with value of β, seems like the code above works good for both zero and non-zero degrees. Also, I didn't care about function domains, you should reproduce the derivation to ensure that all parameters and solution are well-defined.

Wrong solution to a quadratic equation

1 Answers