[SOLVED] Wrong solution to a quadratic equation

Issue

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.

Solution

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.

Answered By – Stepan Zakharov

Answer Checked By – Pedro (BugsFixing Volunteer)