# Fortran Getting the cubic root

1. Sep 23, 2016

### Amany Gouda

I tried to calculate the cubic root by using the method that are exist in Numerical receipes 77 but I got no answer and I don't know my mistake .
Also, I tried by using Cardino method but Also I couldn't success to get an answer.
Can any read my codes and tell me where is my errors or provide me with more ideas.

first code:
Code (Fortran):

! this program is made by using numerical receipes book version 77.

2 !------------------------------------------------------------------

3 program main

4 implicit none

5 real*8::Q,R,Theta,b,c,d

6 real*8,parameter::pi=3.14

7 complex*8::z1,z2,z3,x1,x2,x3

8 print*,'b,c,d'

10

11 Q=((b**2)-(3*c))/9

12 R=((2*(b**3))-(9*b*c)+(27*d))/54

13

14

15 Theta=acos(R/(sqrt(Q**3)))

16

17 x1=-2*(sqrt(Q))*(cos(Theta/3))-(b/3)

18 x2=-2*(sqrt(Q))*(cos((Theta+2*pi)/3))-(b/3)

19 x3=-2*(sqrt(Q))*(cos((Theta-2*pi)/3))-(b/3)

20

21

22 print*,'x1=',x1,'x2=',x2,'x3=',x3

23

24

25 z1= sqrt(x1)

26 z2= sqrt(x2)

27 z3=  sqrt(x3)

28

29

30 print*,' z1,z2,z3',z1,z2,z3

31 end program

The second code:
Code (Fortran):

! this program is made to used for calculating the cardino method for the our system parameters.

2 !-----------------------------------------------------------------------------------------------

3 program main

4 Implicit none

5 Integer::i,a

6 real*8,dimension (100):: Il,Il2,a0,omega,nominator,dominator

7 real*8,dimension(100)::f,lamda,kk,sum0,sum1,b,c,d,cardino3,Delta3

8 real*8,parameter::mass_ratio=0.00870

9 complex*8,dimension(100)::z1,z2,z3

10

11 a=1

12 Il=0.0

13 Il2=0.0

14 sum0=0.0

15 sum1=0.0

16

17 do i =1,100

18

19    Il(i)=sum0(i)+0.1

20

21    omega(i)=2.359*Il(i)

22    nominator=2-omega(i)

23    dominator=1-omega(i)

24    f(i)=nominator(i)/dominator(i)

25    lamda(i)=sqrt(f(i))

26    kk(i)=1/lamda(i)

27

28   Il2(i)=sum1(i)+0.001

29   a0(i)=0.68*Il2(i)

30

31    sum0=Il(i)

32    sum1=Il2(i)

33

34

35   print*,'Il=',Il(i)

36   print*,'Tl2=',Il2(i)

37   print*,'omega=',omega(i)

38   print*,'a0=',a0(i)

39   print*,'lamda=',lamda(i)

40   print*,'kk=',kk(i)

41

42 call eqn_parameter(omega(i),a0(i),kk(i),b(i),c(i),d(i))

43

44    print*,'b=',b(i)
45    print*,'c=',c(i)

46    print*,'d=',d(i)

47

48 call cardino_method(a,b(i),c(i),d(i),cardino3(i),Delta3(i))

49

50   print*,' carino_solution',cardino3(i)

51   print*,'Delta0=',Delta3(i)

52

53 call fsub(cardino3(i),b(i),Delta3(i),z1(i),z2(i),z3(i))

54

55

56     print *, "Omega1_square= ",z1(i)

57     print *, "Omega2_square= ",z2(i)

58     print *, "Omega3_square= ",z3(i)

59

60

61 end do

62 end program

63

64

65 subroutine eqn_parameter(x,t,s,f1,f2,f3)

66     implicit none

67     real*8,dimension(100),intent(in)::x,t,s

68     real*8,dimension(100),intent(out)::f1,f2,f3

69     real*8,dimension(100)::bb,cc1,cc2,cc3,dd

70

71     bb=-(2.0+x)

72     cc1=2.0*x

73     cc2=2.0*s*t*0.00870

74     cc3=(1.0d0+cc1-cc2)

75     dd=-(x+cc2)

76

77 end subroutine eqn_parameter

78

79

(mentor note: code edited to remove bolding and to use code tags)

Last edited by a moderator: Sep 23, 2016
2. Sep 23, 2016

### Staff: Mentor

That's too vague. Can you give an example of input and output?

While not related to your problem, this is an eyesore: real*8,parameter::pi=3.14
I'm sure you can do better precision-wise.

Also, the equivalent to real*8 is complex*16.

3. Sep 23, 2016

### Amany Gouda

I mean no output for the cubic root please check the attached file.
the output for the omega square and omega root is zeros.

#### Attached Files:

• ###### 23092016.txt
File size:
65.3 KB
Views:
40
4. Sep 23, 2016

### Staff: Mentor

In the file, I see
Code (Text):

b=   6.9273636908778262E-310
c=   6.9532115335400158E-310
d=   6.9532620141218920E-310

What is the range of a real*8? What do you thinks happens in lines 11 and 12 of the first program?

5. Sep 23, 2016

### Amany Gouda

For the range, I never think for a range to my result but I thought real*8 represent double precision that mean 64 bit.

DOUBLE PRECISION
REAL (KIND=8)
REAL*8 8 (64 bits) Double-precision real floating-point values in IEEE T_float format ranging from 2.2250738585072013D-308 to 1.7976931348623158D308. Values between 2.2250738585072008D--308 and 4.94065645841246544D--324 are denormalized 3.

for line 11 and 12;
it is a method from the numerical receipes book to get the cubic root by the following steps;

please check page number 179 in the numerical receipes 77;
http://perso.ens-lyon.fr/christophe.winisdoerffer/INTRO_NUM/NumericalRecipesinF77.pdf

I calculated Q and R to get help in calculating the three roots of the cubic equation.

6. Sep 23, 2016

### Staff: Mentor

My point is that your are calculating the square and the cube of a number that is ~10-310. The answer will necessarily be zero.

7. Sep 23, 2016

### Amany Gouda

Then please check this modification and tell me where I should consider the range in my new modified code. also, why in matlab the same equation produced full results without considering the range.
! this program is made to used for calculating the cardino method for the our system parameters.

2 !-----------------------------------------------------------------------------------------------

3 program main

4 Implicit none

5 Integer::i,a
6 real*8,dimension (100):: Il,Il2,a0,omega,nominator,dominator

7 real*8,dimension(100)::f,lamda,kk,sum0,sum1,b,c,d,constant_c1,Delta3

8 real*8,parameter::mass_ratio=0.00870

9 complex*8,dimension(100)::z1,z2,z3

10 complex*8,dimension(100)::zz1,zz2,zz3

11

12

13 a=1

14 Il=0.0

15 Il2=0.0

16 sum0=0.0

17 sum1=0.0

18

19 do i =1,100

20

21 Il(i)=sum0(i)+0.1

22

23 omega(i)=2.359*Il(i)

24 nominator=2-omega(i)

25 dominator=1-omega(i)

26 f(i)=nominator(i)/dominator(i)

27 lamda(i)=sqrt(f(i))

28 kk(i)=1/lamda(i)

29

30 Il2(i)=sum1(i)+0.001

31 a0(i)=0.68*Il2(i)

32

33 sum0=Il(i)

34 sum1=Il2(i)

35

36

37 print*,'Il=',Il(i)

38 print*,'Tl2=',Il2(i)

39 print*,'omega=',omega(i)

40 print*,'a0=',a0(i)

41 print*,'lamda=',lamda(i)

42 print*,'kk=',kk(i)
44 call eqn_parameter(omega(i),a0(i),kk(i),b(i),c(i),d(i))

45

46 print*,'b=',b(i)

47 print*,'c=',c(i)

48 print*,'d=',d(i)

49

50 call constantc(a,b(i),c(i),d(i),constant_c1(i),Delta3(i))

51

52 print*,' constant_c',constant_c1(i)

53 print*,'Delta0=',Delta3(i)

54

55 call fsub(constant_c1(i),b(i),Delta3(i),z1(i),z2(i),z3(i))

56

57

58 print *, "Omega1_square= ",z1(i)

59 print *, "Omega2_square= ",z2(i)

60 print *, "Omega3_square= ",z3(i)

61

62

63 zz1(i)=sqrt(z1(i))

64 zz2(i)=sqrt(z2(i))

65 zz3(i)=sqrt(z3(i))

66

67 print *, "Omega1_root= ",zz1(i)

68 print *, "Omega2_root= ",zz2(i)

69 print *, "Omega3_root= ",zz3(i)

70

71 end do

72

73 end program

74

75

76 subroutine eqn_parameter(x,t,s,f1,f2,f3)

77 implicit none

78 real*8,dimension(100),intent(in)::x,t,s

79 real*8,dimension(100),intent(out)::f1,f2,f3

80 real*8,dimension(100)::bb,cc1,cc2,cc3,dd

81

82 bb=-(2.0+x)

83 cc1=2.0*x

84 cc2=2.0*s*t*0.00870

85 cc3=(1.0d0+cc1-cc2)

86 dd=-(x+cc2)
85 cc3=(1.0d0+cc1-cc2)

86 dd=-(x+cc2)

87

88 end subroutine eqn_parameter

89

90

91 subroutine constantc(aa,be,ci,di,constant_c,Delta0)

92

93 implicit none

94 real*8,dimension(100),intent(in)::be,ci,di

95 real*8,dimension (100),intent(out)::constant_c,Delta0

96 integer::aa

97 real*8,dimension(100)::Delta1,m,BB,mm,CC,NN,cardino1

98 double precision::THIRD= 0.3333333333333333333333333333333333333333333333333333333333333333333333333333333333D0

99

100

101 Delta0= (be**2)-(3*aa*ci)

102 Delta1= (2*(be**3))-(9*aa*be*ci)+(27*(aa**2)*di)

103 m=(Delta0)**3

104 BB=(Delta1)**2

105 mm= BB-(4*m)

106 CC=sqrt(mm)

107 NN=(Delta1+CC)/2

108 constant_c= SIGN((ABS(NN))**THIRD, NN)

109

110

111 end subroutine constantc

112

113 subroutine fsub(g,w,q,u,tt,r)

114

115 implicit none

116 real*8,intent(in),dimension(100)::g,w,q

117 complex*8,dimension(100),intent(out)::u,tt,r

118 complex*8,dimension(100)::z1,z2,z3

119

120 z1=(1,0)

121 u=(-1/3)*(q+(z1*g)+(w/(z1*g)))

122 z2=(-0.5,0.866)

123 tt=(-1/3)*(q+(z2*g)+(w/(z2*g)))

124 z3=(-0.5,-0.866)

125 r=(-1/3)*(q+(z3*g)+(w/(z3*g)))

126 end subroutine fsub

127