Calculating Cubic Root: Numerical Recipes & Cardino Method

  • Context: Fortran 
  • Thread starter Thread starter Amany Gouda
  • Start date Start date
  • Tags Tags
    Cubic Root
Click For Summary
SUMMARY

The forum discussion centers on the challenges faced in calculating cubic roots using the Numerical Recipes 77 and Cardano methods. The user reported issues with obtaining valid outputs from their Fortran code implementations, specifically noting that the results for cubic roots were returning as zero. Key points include the importance of understanding the range of real*8 data types in Fortran, which represent double-precision floating-point values, and the necessity of ensuring that input values are within a computable range to avoid underflow issues.

PREREQUISITES
  • Understanding of Fortran programming, specifically version 77.
  • Familiarity with numerical methods for solving polynomial equations, particularly the Cardano method.
  • Knowledge of complex numbers and their representation in Fortran.
  • Comprehension of floating-point precision and the implications of using real*8 data types.
NEXT STEPS
  • Review the Numerical Recipes 77 documentation, focusing on cubic root calculations.
  • Learn about the implications of floating-point precision and how to avoid underflow in numerical computations.
  • Explore advanced debugging techniques in Fortran to identify and resolve issues in code execution.
  • Investigate alternative numerical libraries or tools, such as MATLAB, for cubic root calculations and compare results.
USEFUL FOR

This discussion is beneficial for software developers, numerical analysts, and researchers working with polynomial equations, particularly those using Fortran for scientific computing and requiring accurate root calculations.

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

  • · Replies 7 ·
Replies
7
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 14 ·
Replies
14
Views
4K
  • · Replies 6 ·
Replies
6
Views
5K
  • · Replies 17 ·
Replies
17
Views
5K
  • · Replies 1 ·
Replies
1
Views
6K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 4 ·
Replies
4
Views
2K
  • · Replies 11 ·
Replies
11
Views
7K
  • · Replies 1 ·
Replies
1
Views
2K