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

Hi,

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?

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:

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

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.