Numerical resolution using cardan method

1

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…]

2

Hi,

This has been discussed here.

1 Like
3

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.

2 Likes
4

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).

5

Which Fortran compiler and version are you using?

I’ve slightly adapted the code by @ELNS and @stavros to use a subroutine instead of a derived type:

!
! 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

Example:

$ 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:

1 Like
6

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.

7

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. Algorithm 954: An Accurate and Efficient Cubic and Quartic Equation Solver for Physical Applications: ACM Transactions on Mathematical Software: Vol 41, No 4

A Fortran 90 module is available in the supplemental material.

8

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.