Calculating Cubic Root: Numerical Recipes & Cardino Method

  • Context: Fortran 
  • Thread starter Thread starter Amany Gouda
  • Start date Start date
  • Tags Tags
    Cubic Root
Amany Gouda
Messages
29
Reaction score
0
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:
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'

  9 read*,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:
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:
on Phys.org
Amany Gouda said:
I got no answer
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.
 
DrClaude said:
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.
I mean no output for the cubic root please check the attached file.
the output for the omega square and omega root is zeros.
 

Attachments

In the file, I see
Code:
 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?
 
DrClaude said:
In the file, I see
Code:
 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?
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.
 
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.
 
DrClaude said:
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.

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