Fortran Calculating Cubic Root: Numerical Recipes & Cardino Method

  • Thread starter Thread starter Amany Gouda
  • Start date Start date
  • Tags Tags
    Cubic Root
AI Thread Summary
The discussion centers on difficulties in calculating cubic roots using the Numerical Recipes and Cardano methods, with the user reporting no output from their code. Key issues identified include the use of very small values leading to underflow, resulting in outputs close to zero for coefficients b, c, and d. Participants suggest checking the input values and the calculations in lines 11 and 12 of the first program, as they may be producing negligible results. The user is encouraged to consider the range of real*8 variables and the implications of using such small numbers in calculations. The conversation highlights the need for careful handling of numerical precision in programming for accurate results.
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:
Technology news 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
 

Similar threads

Back
Top