Suppose the following transcendental equation $$ F(u)=\sin u(A-4u^{2})+u(B+Cu)-D=0, $$ where $A=\pi^{2}(x^{2}+y^{2})$, $B=2\pi(1-x^{2}-y^{2})$, $C=4y$, $D=\pi^{2}y$, are the constant values, and $x\in[-1,1${]}, $y\in[-1,1]$. $F(u)$ has multiple roots ($m\leq4$) on interval $u\in[-\pi,\pi]$.
x = 0.5
y = 0.5
u = -pi:0.01:pi
A = pi^2 * (x^2 + y^2)
B = 2 * pi * (1 - x^2 - y^2)
C = 4.0 * y
D = pi^2*y;
F = sin(u).*(A - 4*u.^2) + u.*(B + u*C) - D
plot(u, F)
In my opinion there is no analytic solution of this equation.
Questions:
Which method for solving transcendental equations with multiple roots do you recommend?
Is there any method avoiding the complex number arithmetic?
Solution including complex numbers:
Using the solution described in paper, where $R=\pi$, $h=0$, the following equations
$$ \left(\begin{array}{cccc} G_{1} & G_{2} & G_{3} & G_{4}\\ G_{2} & G_{3} & G_{4} & G_{5}\\ G_{3} & G_{4} & G_{5} & G_{6}\\ G_{4} & G_{5} & G_{6} & G_{7} \end{array}\right)\left(\begin{array}{c} c_{0}\\ c_{1}\\ c_{2}\\ c_{3} \end{array}\right)=-\left(\begin{array}{c} G_{5}\\ G_{6}\\ G_{7}\\ G_{8} \end{array}\right), $$
are solved
$$ c=\left(\begin{array}{c} -0.0371595\\ 0.228188\\ -0.101362\\ -0.912752 \end{array}\right), $$
and
$$ F(u)=-0.0371595+0.228188u-0.101362u^{2}-0.912752u^{3}+u^{4}, $$ with roots $u_{1}=-1.570800$, $u_{2}=-0.666522$, $u_{3}=1.570800$, $u_{4}=2.20097$. This method works well, however, due to the non-trivial implementation of integrals, complex numbers, etc., it is not easy to code in C++.
Mathematica solution:
x = 0.5
y = 0.5
R = Pi;
h = 0;
(*Substitution*)
z = Exp[I*theta];
u = R * z + h;
a = Pi^2*(x^2 + y^2)
b = 2*Pi*(1 - x^2 - y^2)
c = 4.0*y
d = Pi^2*y;
F = Sin[u]*(a - 4*u^2) + u*(b + u*c) - d;
G = 1/F;
(*Integration*)
k = 1;
G1 = N[NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}], 8];
k = 2;
G2 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
k = 3;
G3 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
k = 4;
G4 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
k = 5;
G5 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
k = 6;
G6 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
k = 7;
G7 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
k = 8;
G8 = NIntegrate[G*Exp[I*k*theta], {theta, 0.00, 2*Pi}];
A = {{Re[G1], Re[G2], Re[G3], Re[G4]}, {Re[G2], Re[G3], Re[G4],
Re[G5]}, {Re[G3], Re[G4], Re[G5], Re[G6]}, {Re[G4], Re[G5],
Re[G6], Re[G7]}};
h = -{{Re[G5]}, {Re[G6]}, {Re[G7]}, {Re[G8]}};
c = Inverse[A].h
(*Solve cubic equation*)
f = c[[1, 1]] + c[[2, 1]]*X + c[[3, 1]]*X^2 + c[[4, 1]]*X^3 + X^4;
S = Solve[f == 0, X];
s = S[[All, 1, 2]];
(*Find roots*)
u1 = s[[1]] * R
u2 = s[[2]] * R
u3 = s[[3]] * R
u4 = s[[4]] * R