I need to solve an equation of the type
a(x)*x**4 + b(x)**2 + c(x) = 0
where a(x), b(x), and c(x) are functions of the type
f(x) = c1 / (c2 + c3*x)
I’m looking for some Fortran routine/function/module to help me.
I found many powerful codes to solve systems of equations, but using them for a single equation feels like using a cannon to kill a fly.
Could you recommend any option?
Thanks,
Could Newton’s method be of help here?
Hi @mcb, welcome to Discourse.
One option would be this package:
Looking at the following two resources:
I couldn’t find any special solvers for rational functions.
If I’m not mistaken, you can rewrite your equation,
a_1/(a_2 + a_3 x) x^4 + (b_1/(b_2 + b_3 x))^2 + c_1/(c_2 + c_3 x)
like this
f(x) = P(x) / Q(x)
P(x) = (a_1 x^4 (b_3 x + b_2)^2 (c_3 x + c_2) + (a_3 x + a_2) (b_1^2 (c_3 x + c_2) + c_1 (b_3 x + b_2)^2))
Q(x) = ((a_3 x + a_2) (b_3 x + b_2)^2 (c_3 x + c_2))
Now the roots of Q(x) are the poles, and the roots of P(x) are the roots. Is that correct? What happens when the roots and poles are at the same x?
For polynomial functions special solvers exist. I would recommend polyroots:
Following on from @ivanpribec , you can expand your equation as
a2*b2**2*c1 + a2*b1**2*c2 + (a3*b2**2*c1 + 2*a2*b2*b3*c1 + a3*b1**2*c2 + a2*b1**2*c3)*x+(2*a3*b2*b3*c1 + a2*b3**2*c1 + a3*b1**2*c3)*x**2 + a3*b3**2*c1*x**3 + a1*b2**2*c2*x**4 + (2*a1*b2*b3*c2 + a1*b2**2*c3)*x**5 + (a1*b3**2*c2 + 2*a1*b2*b3*c3)*x**6 + a1*b3**2*c3*x**7
which is a 7’th order polynomial.
Then construct the companion matrix to this polynomial. Pass this matrix into the LAPACK routine dgeev to compute the (possibly complex) eigenvalues. These are all the roots of your original equation.
Minus the ones that potentially align with the poles.
I’ve posted an example of using the companion matrix before for cubic polynomials: Cardano's solution of the cubic equation - #5 by ivanpribec
Thanks for the input.
If I had to code it myself, Newton’s method would be my first attempt.
However, I forgot to mention that this equation must be solved thousands of times in each program execution.
Therefore, I’m seeking a professional solution with the fastest algorithms and best i/o and memory management. Developing that at this optimal level goes beyond my expertise. And I’m sure that the Fortran community should already have excellent solutions.
Thanks, @ivanpribec and @jkd2022, for your comments.
It’s embarrassing I didn’t realize that I could rewrite my equation as a high-order polynomial! Do you know when you’re so deep into the box that you miss entirely the outside perspective?
Indeed, things should be much simpler for a polynomial equation. I will check the packages you mentioned.
Concerning poles, although the program should detect them occasionally, my specific problem is defined in a domain where the solutions should be in well-behaved regions.
Thanks again,
Does it need to solved thousands of time concurrently or sequentially?
In case it’s the former, you could try a vectorized approach. Newton’s method is suitable for this, assuming your problem has some local convergence guarantees. The SciPy optimize.newton routine supports this:
import numpy.array_api as xp
from scipy.optimize import newton
def r(x,c):
return c[0] / (c[1] + c[2]*x)
def rp(x,c):
return -c[0]*c[2] / (c[1] + c[2]*x)**2
def f(x,a,b,c):
return r(x,a) * x**4 + (r(x,b))**2 + r(x,c)
# d/dx(a(x) x^4 + b(x)^2 + c(x)) = x^4 a'(x) + 4 x^3 a(x) + 2 b(x) b'(x) + c'(x)
def fp(x):
return x**3 * (rprime(x,a) + 4 * r(x,a)) + 2* r(x,b) * rprime(x,b) + rprime(x,c)
if __name__ == "__main__":
a = ( a1, a2, a3 )
b = ( b1, b2, b3 )
c = ( c1, c2, c3 )
x0 = xp.random.rand(1000)
root = newton(f,x0,fprime,args=(a,b,c))
It would take a some effort to replicate this in Fortran.
What counts as fast for you? On a modern (Intel) CPU the classic zeroin routine from the 1970’s can solve a 10^5 root-solving problems in 0.02 s (https://twitter.com/IvanPribec/status/1741970602435526885). That is 0.2 µs per solve. This is for an extremely simple function, so we can assume the time is mostly spent in the solver.
Maybe you’d like to give ForSolver a try.
Here’s an example using the Newton method:
program newton
use kinds
use forsolver
implicit none
type(nlsolver) :: nls
real(rk) :: x
call nls%set_options(nl_method='newton', maxit=100, TolFun=1e-6_rk)
call nls%solve(F=F, dFdx=dFdx, x0=0.0_rk, x_sol=x)
contains
function r(x,c1,c2,c3)
real(rk), intent(in) :: x
real(rk), intent(in) :: c1, c2, c3
real(rk) :: r
r = c1 / (c2 + c3*x)
end function
function drdx(x,c1,c2,c3)
real(rk), intent(in) :: x
real(rk), intent(in) :: c1, c2, c3
real(rk) :: drdx
drdx = -c1 * c3 / (c2 + c3 * x)**2
end function
function F(x)
real(rk), intent(in) :: x
real(rk) :: F
real(rk) :: c1, c2, c3
c1 = 10.0_rk
c2 = 12.0_rk
c3 = 30.0_rk
F = r(x,c1,c2,c3) * x**4 + r(x,c1,c2,c3)**2 + r(x,c1,c2,c3)
end function
function dFdx(x)
real(rk), intent(in) :: x
real(rk) :: c1, c2, c3
real(rk) :: dFdx
real(rk) :: dr_dx
c1 = 10.0_rk
c2 = 12.0_rk
c3 = 30.0_rk
dFdx = drdx(x,c1,c2,c3) * x**4 + r(x, c1, c2, c3) * 4.0_rk * x**3&
+ 2.0_rk*drdx(x,c1,c2,c3)*r(x, c1, c2, c3)&
+ drdx(x,c1,c2,c3)
end function
end program
-----------------------------------------------
maxit x0 tol
100 0.00000000 0.1000E-05
-----------------------------------------------
start newton
-----------------------------------------------
it xn F(xn) dF(xn)/dxn
0 0.0000 0.1528E+01 -0.5556E+01
1 0.2750 0.7405E+00 -0.1417E+01
2 0.7975 0.4684E+00 0.1089E+00
3 -3.5056 -0.1631E+02 0.1325E+02
4 -2.2749 -0.4908E+01 0.5772E+01
5 -1.4246 -0.1559E+01 0.2344E+01
6 -0.7592 -0.3751E+00 0.2978E+01
7 -0.6332 0.3838E+00 0.1186E+02
8 -0.6656 0.7384E-01 0.7689E+01
9 -0.6752 0.4195E-02 0.6838E+01
10 -0.6758 0.1534E-04 0.6788E+01
11 -0.6758 0.2069E-09 0.6788E+01
-----------------------------------------------
end newton
-----------------------------------------------
x_sol = -0.67580449156891309
And here’s the same example using the complex step Newton method:
program newton_quasi_cs
use kinds
use forsolver
implicit none
type(nlsolver) :: nls
complex(rk) :: x
call nls%set_options(nl_method='newton-quasi-cs', maxit=100, TolFun=1e-6_rk, cs_tol=tiny(0.0_rk))
call nls%solve(F=F, x0=(0.0_rk, 0.0_rk), x_sol=x)
contains
function r(x,c1,c2,c3)
complex(rk), intent(in) :: x
complex(rk), intent(in) :: c1, c2, c3
complex(rk) :: r
r = c1 / (c2 + c3*x)
end function
function F(x)
complex(rk), intent(in) :: x
complex(rk) :: F
complex(rk) :: c1, c2, c3
c1 = (10.0_rk, 0.0_rk)
c2 = (12.0_rk, 0.0_rk)
c3 = (30.0_rk, 0.0_rk)
F = r(x,c1,c2,c3) * x**4 + r(x,c1,c2,c3)**2 + r(x,c1,c2,c3)
end function
end program
-----------------------------------------------
maxit x0 tol
100 0.00000000 0.1000E-05
-----------------------------------------------
start newton
-----------------------------------------------
it xn F(xn) dF(xn)/dxn
0 0.0000 0.1528E+01 -0.5556E+01
1 0.2750 0.7405E+00 -0.1417E+01
2 0.7975 0.4684E+00 0.1089E+00
3 -3.5056 -0.1631E+02 0.1325E+02
4 -2.2749 -0.4908E+01 0.5772E+01
5 -1.4246 -0.1559E+01 0.2344E+01
6 -0.7592 -0.3751E+00 0.2978E+01
7 -0.6332 0.3838E+00 0.1186E+02
8 -0.6656 0.7384E-01 0.7689E+01
9 -0.6752 0.4195E-02 0.6838E+01
10 -0.6758 0.1534E-04 0.6788E+01
11 -0.6758 0.2069E-09 0.6788E+01
-----------------------------------------------
end newton
-----------------------------------------------
x_sol = -0.67580449156891309
It looks like you have already gotten some good replies. I just wanted to add that sometimes the best algorithm choice depends on what is known about the problem. Are all the solutions of interest known to be purely real, or purely imaginary? Sometimes, only the solution(s) within a predefined interval are wanted, so any complications (such as poles, or multiple roots) outside of that interval can be eliminated. Sometimes, there are rigorous lower and upper bounds that are known – this suggests a combination of methods might be useful, such as a combination of Newton and bisection (e.g. bisection when the Newton steps are outside the range, and Newton steps when within the neighborhood of second-order convergence). If rigorous bounds are not available, then approximate bounds might be generated using trust-radius methods. If you are solving for many solutions, can information from one be used to advantage in the next – this could include the computation of shared intermediates, or starting guesses for the iterative procedures.
In the triad of “better, faster, cheaper”, which is/are the most important? One can often get two of those at the expense of the third.