Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran Getting the cubic root

  1. Sep 23, 2016 #1
    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'

      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:
    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. jcsd
  3. Sep 23, 2016 #2

    DrClaude

    User Avatar

    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.
     
  4. Sep 23, 2016 #3

    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:

  5. Sep 23, 2016 #4

    DrClaude

    User Avatar

    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?
     
  6. Sep 23, 2016 #5

    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.
     
  7. Sep 23, 2016 #6

    DrClaude

    User Avatar

    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.
     
  8. Sep 23, 2016 #7
    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
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Getting the cubic root
  1. Cubic Root (Replies: 3)

  2. Root Finding (Replies: 5)

  3. Root programming (Replies: 3)

  4. Root programming (Replies: 3)

Loading...