How to resolve the algorithm Roots of a quadratic function step by step in the Fortran programming language

Published on 12 May 2024 09:40 PM

How to resolve the algorithm Roots of a quadratic function step by step in the Fortran programming language

Table of Contents

Problem Statement

Write a program to find the roots of a quadratic equation, i.e., solve the equation

a

x

2

b x + c

0

{\displaystyle ax^{2}+bx+c=0}

. Your program must correctly handle non-real roots, but it need not check that

a ≠ 0

{\displaystyle a\neq 0}

. The problem of solving a quadratic equation is a good example of how dangerous it can be to ignore the peculiarities of floating-point arithmetic. The obvious way to implement the quadratic formula suffers catastrophic loss of accuracy when one of the roots to be found is much closer to 0 than the other. In their classic textbook on numeric methods Computer Methods for Mathematical Computations, George Forsythe, Michael Malcolm, and Cleve Moler suggest trying the naive algorithm with

a

1

{\displaystyle a=1}

,

b

10

5

{\displaystyle b=-10^{5}}

, and

c

1

{\displaystyle c=1}

. (For double-precision floats, set

b

10

9

{\displaystyle b=-10^{9}}

.) Consider the following implementation in Ada: As we can see, the second root has lost all significant figures. The right answer is that X2 is about

10

− 6

{\displaystyle 10^{-6}}

. The naive method is numerically unstable. Suggested by Middlebrook (D-OA), a better numerical method: to define two parameters

q

a c

/

b

{\displaystyle q={\sqrt {ac}}/b}

and

f

1

/

2 +

1 − 4

q

2

/

2

{\displaystyle f=1/2+{\sqrt {1-4q^{2}}}/2}

and the two roots of the quardratic are:

− b

a

f

{\displaystyle {\frac {-b}{a}}f}

and

− c

b f

{\displaystyle {\frac {-c}{bf}}}

Task: do it better. This means that given

a

1

{\displaystyle a=1}

,

b

10

9

{\displaystyle b=-10^{9}}

, and

c

1

{\displaystyle c=1}

, both of the roots your program returns should be greater than

10

− 11

{\displaystyle 10^{-11}}

. Or, if your language can't do floating-point arithmetic any more precisely than single precision, your program should be able to handle

b

10

6

{\displaystyle b=-10^{6}}

. Either way, show what your program gives as the roots of the quadratic in question. See page 9 of "What Every Scientist Should Know About Floating-Point Arithmetic" for a possible algorithm.

Let's start with the solution:

Step by Step solution about How to resolve the algorithm Roots of a quadratic function step by step in the Fortran programming language

Source code in the fortran programming language

PROGRAM QUADRATIC

 IMPLICIT NONE
 INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15)
 REAL(dp) :: a, b, c, e, discriminant, rroot1, rroot2
 COMPLEX(dp) :: croot1, croot2

 WRITE(*,*) "Enter the coefficients of the equation ax^2 + bx + c"
 WRITE(*, "(A)", ADVANCE="NO") "a = "
 READ *, a
 WRITE(*,"(A)", ADVANCE="NO") "b = "
 READ *, b
 WRITE(*,"(A)", ADVANCE="NO") "c = "
 READ *, c
 
 WRITE(*,"(3(A,E23.15))") "Coefficients are: a = ", a, "   b = ", b, "   c = ", c
 e = 1.0e-9_dp
 discriminant = b*b - 4.0_dp*a*c
 
 IF (ABS(discriminant) < e) THEN
    rroot1 = -b / (2.0_dp * a)
    WRITE(*,*) "The roots are real and equal:"
    WRITE(*,"(A,E23.15)") "Root = ", rroot1
 ELSE IF (discriminant > 0) THEN
    rroot1 = -(b + SIGN(SQRT(discriminant), b)) / (2.0_dp * a)
    rroot2 = c / (a * rroot1)
    WRITE(*,*) "The roots are real:"
    WRITE(*,"(2(A,E23.15))") "Root1 = ", rroot1, "  Root2 = ", rroot2
 ELSE
    croot1 = (-b + SQRT(CMPLX(discriminant))) / (2.0_dp*a) 
    croot2 = CONJG(croot1)
    WRITE(*,*) "The roots are complex:" 
    WRITE(*,"(2(A,2E23.15,A))") "Root1 = ", croot1, "j ", "  Root2 = ", croot2, "j"
 END IF


COMPUTE ROOTS OF A QUADRATIC FUNCTION - 1956
      READ 100,A,B,C
 100  FORMAT(3F8.3)
      PRINT 100,A,B,C
      DISC=B**2-4.*A*C
      IF(DISC),1,2,3
 1    XR=-B/(2.*A)
      XI=SQRT(-DISC)/(2.*A)
      XJ=-XI
      PRINT 311
      PRINT 312,XR,XI,XR,XJ
 311  FORMAT(13HCOMPLEX ROOTS)
 312  FORMAT(4HX1=(,2E12.4,6H),X2=(,2E12.4,1H))
      GO TO 999
 2    X1=-B/(2.*A)
      X2=X1
      PRINT 321
      PRINT 332,X1,X2
 321  FORMAT(16HEQUAL REAL ROOTS)
      GO TO 999
 3    X1= (-B+SQRT(DISC)) / (2.*A)
      X2= (-B-SQRT(DISC)) / (2.*A)
      PRINT 331
      PRINT 332,X1,X2
 331  FORMAT(10HREAL ROOTS)
 332  FORMAT(3HX1=,E12.5,4H,X2=,E12.5)
 999  STOP


  

You may also check:How to resolve the algorithm Multiple regression step by step in the FreeBASIC programming language
You may also check:How to resolve the algorithm Literals/String step by step in the Picat programming language
You may also check:How to resolve the algorithm Keyboard input/Obtain a Y or N response step by step in the Inform 7 programming language
You may also check:How to resolve the algorithm Price fraction step by step in the OCaml programming language
You may also check:How to resolve the algorithm Bulls and cows/Player step by step in the zkl programming language