Good evening dear all, I need your expertise. Please help me find a Fortran code that solves a third degree polynomial, with real coefficients a, b, c, d and a real parameter m to vary. Thanks verry much for your response. [Uploading: 16026184485728987933868987474674.jpg…]
This has been discussed here.
Hello Justin, welcome to Discourse.
If you are asking for a solution of the equation
am^3 + bm^2 + cm + d=0
the answer is indeed in the link posted by @pcosta.
Thank you very much, I am talking about an equation of the following form : am^3 + (b+k)m^2 + (k-c)m + d = 0, where k € [1, 4] for example. Also I tried using the code published previously for the simple case but it shows also error (solve has no implicit type).
Which Fortran compiler and version are you using?
! ! Cardano's method for cubic polynomials ! ! Adapted from the example by EverLookNeverSee: ! https://github.com/EverLookNeverSee/FCS/blob/master/src/numerical%20methods/Exercise_14.f90 ! ! Related discussion is found at: ! https://fortran-lang.discourse.group/t/cardanos-solution-of-the-cubic-equation/111/6 ! module m_cardano implicit none private public :: wp, solve integer, parameter :: wp = kind(1.0e0) contains subroutine solve(a, b, c, d, r1, r2, r3) ! declaring dummy parameters and local variables real(wp), intent(in) :: a, b, c, d real(wp), intent(out) :: r1 complex(wp), intent(out) :: r2, r3 real(wp) :: Q, R, S, T, r_part, i_part, temp ! calculating values of Q and R Q = (3.0_wp*a*c - b**2) / (9.0_wp * a**2) R = (9.0_wp*a*b*c - 27.0_wp*d* a**2 - 2.0_wp * b**3) / (54.0_wp * a**3) ! using temp adjunct variable and abs intrinsic function in order to ! prevent `negative real to a real power` error temp = R + sqrt(Q**3 + R**2) S = sign(1.0_wp, temp) * abs(temp)**(1.0_wp / 3.0_wp) temp = R - sqrt(Q**3 + R**2) T = sign(1.0_wp, temp) * abs(temp)**(1.0_wp / 3.0_wp) ! calculating real and imaginary parts of complex roots r_part = -(s + t) / 2.0_wp - b / (3.0_wp * a) i_part = (sqrt(3.0_wp) / 2.0_wp) * (s - t) ! calculating and returning real and complex roots r1 = S + T - b / (3.0_wp * a) r2 = cmplx(r_part, i_part) r3 = cmplx(r_part, -i_part) end subroutine end module program test_cardano use m_cardano implicit none real(wp) :: a, b, c, d real(wp) :: r1 complex(wp) :: r2, r3 print *, "Cardano's solution for solving cubic equation" print *, "Cubic equation: ax^3 + bx^2 + cx + d = 0" print *, "Enter coefficients: a, b, c, d respectively:" ! getting coefficients of cubic equation using user input read *, a, b, c, d ! calling solve function call solve(a,b,c,d,r1,r2,r3) ! printing results print *, "Real root:", r1 print *, "Complex roots:", r2, r3 end program
$ gfortran -Wall m_cardano.f90 $ ./a.out Cardano's solution for solving cubic equation Cubic equation: ax^3 + bx^2 + cx + d = 0 Enter coefficients: a, b, c, d respectively: 4. 3. 2. 1. Real root: -0.605829597 Complex roots: (-7.208520174E-02,0.638326645) (-7.208520174E-02,-0.638326645)
Using the following SymPy commands I can verify the results:
$ python >>> import sympy as sp >>> x = sp.symbols('x') >>> eq = 4*x**3 + 3*x**2 + 2*x + 1 >>> res = sp.roots(eq) >>> resn = [key.evalf() for key in res.keys()] [-0.605829586188268, -0.072085206905866 - 0.638326735148376*I, -0.072085206905866 + 0.638326735148376*I]
The tiny difference from the Fortran results is likely due to the use of single precision in the Fortran example.
Maybe it is worth pointing out, the implementation given above is sensitive to round-off errors. A quick search for a robust cubic equation solver gave the following pages which might be helpful:
A numerically robust algorithm for cubic roots can also be found in the work:
Herbison-Evans, D. (1995). Solving quartics and cubics for graphics. In Graphics Gems V (pp. 3-15). Academic Press.
A recent TOMS paper, on the other hand, uses an optimized Newton-Raphson procedure:
Flocke, N. (2015). Algorithm 954: An accurate and efficient cubic and quartic equation solver for physical applications. ACM Transactions on Mathematical Software (TOMS) , 41 (4), 1-24. https://doi.org/10.1145/2699468
A Fortran 90 module is available in the supplemental material.
I am using Fortran 95 on Ubuntu, I will use the different codes and the different basic documents to which you have directed me. It is really nice. Thank you very much for your generosity.