How to resolve the algorithm Roots of a quadratic function step by step in the Fortran programming language
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