# CHALLENGE using root\_fortran to find root of a complex function

**URL:** <https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489>\
**Category:** Uncategorized\
**Created:** [September 8, 2023, 3:07pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489 "2023-09-08T15:07:32Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Mubaraqlanre](https://avatars.discourse-cdn.com/v4/letter/m/7bcc69/32.png) [@Mubaraqlanre](https://fortran-lang.discourse.group/u/Mubaraqlanre)\
**Post date:** [September 8, 2023, 3:07pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/1 "2023-09-08T15:07:32Z")

</div>

Greetings everyone

I am working on a modeling project that involves some complex computations. I transitioned a model from R to Fortran due to performance concerns that hindered the efficiency of the original model being to slow in R. Within this Fortran environment, I am in need of robust root-finding algorithms. While exploring, I came across the use of root-fortran library for finding root in modern Fortran. However on translating the r code to fortran and solving the code using similar tolerance as used in R code I have not been getting the same values as in the R code. The R function is a bit complex and translated to fortran as below

> function Conc\_calc(zHV) result(result)  
> real(Kind=8), intent(in) :: zHV  
> real(Kind=8) :: A, B, Kw, K1, K2, K3, H1, H2, zH  
> real(Kind=8) :: zH0, pH0, pHI, buff\_L  
> real(kind=8), dimension(12) :: vec1, Vec2, vec3, vec4, vec5, vec6, vec7  
> real(kind=8), dimension(12) :: vec8, vec9, vec10, vec11, vec12  
> real(kind=8) :: Amatrix(12,12), Bvector(12)  
> real(kind=8) :: roots(12)  
> real(kind=8) :: result
> 
> ```
> integer :: n=12 ! specify the size of the matrix
>     
>     
> A = x10 + x20 + x30 + InN - OutN
>     
> B = y10 + y20 + y30 + y40 + InC
>     
>                
>     
> Kw = 10.0**(-14.0) ! Kw = 10.0**(-14.0) at 25 C
>     
>     
> !Print *, 'zHV', Tf, ti
>                                     
> K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15))) !defines the first H2CO3 dissociation into HCO3-
>     
> K2 = (10.0**(-1.0 * (((2902.39 / Tf) + (0.02379 * Tf) - 6.4980)))) !defines the second HCO3- dissociation into CO3--
>     
> K3 = (10.0**(-1.0 * (((3404.71 / Tf) + (0.032786 * Tf) - 14.8435)))) !defines the third H2CO3-- dissociation into H2O and co2                   
>     
> H1 = (56.0 * exp(4092.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * (1.0 / Tf) * &
>     
> (1.013 * 10.0 ** 5)) 
>     
> H2 = ((1.0 / 29.41) * exp(2400.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * &
>     
> (1.0 / Tf) * (1.013 * 10.0 ** 5))
> 
>                                       
>     
> zH = zHV * V
> 
> if(ti == 2 .and. ti3 == 3) then
> write(*,*) 'zH is ',zH, ti, ti3                        
> end if  
>     
> vec1 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, -zH, 0.0_wp, & 
> 0.0_wp, (1.0_wp + H1 * (V / Vg)), 0.0_wp/) ! x1
>     
> vec2 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, K1* V, 0.0_wp, 0.0_wp, H1 * (V / Vg), 0.0_wp/) ! x2
> vec6 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
> vec3 = (/0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, -zH, 0.0_wp, 0.0_wp, H2* (V / Vg)/) ! y1
>     
> vec4 = (/0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, K2*V, -zH, 0.0_wp, H2 * (V / Vg)/) ! y2
>     
> vec5 = (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &  
> K3* V, 0.0_wp, (1.0_wp + H2 * (V / Vg))/) ! y3
>     
> vec7 = (/0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! y4
>     
> vec8 = (/1.0_wp, 0.0_wp, 0.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH2
>     
> vec9 = (/-1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK3
>     
> vec10 = (/0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK2
>     
> vec11 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK1
>     
> vec12 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH1
>     
>                 
>     
> Amatrix = reshape([vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, &
>     
> vec11, vec12], shape(Amatrix))
>     
> ! Set the column names
>     
> !character(2) :: colnames(12) = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", &
>     
> ! "RK2", "RK1", "RH1"]
>     
> Bvector = (/y30, y20 + InC, y10, y40 - OutC, x20 + InN, x10, x30 - OutN, 0.0_wp, 0.0_wp, 0.0_wp, &
> H1 * A * (V / Vg), H2 * B * (V / Vg)/) ! const
>     
> call solve_system(n, Amatrix, Bvector, roots) !I called lapack to solve the system         
>     
> zH0 = zHV0 * V
> pH0 = -log10(zH0 / V)
> pHI = -log10(zH / V) !I renamed the line pH to PHI
> result =(zH0 - InC + (-roots(9) + roots(10) + roots(11)) + (buff_L * (pHI - pH0))) - ((zH))
> end function Conc_calc
> 
> ```

I used the root module of root fortran and I called it as below

> call root\_scalar(‘brent’,Conc\_calc, 0.0\_wp, 1.0\_wp, zHV, result, iflag, atol=1.0\_wp,Ftol = 0.10\_wp\*\*(0.25\_wp))

However for the last 3 days I have been stuck as I havent been getting the same result as the R code. I was wondering if I got a step wrong or not so I started to print every single variable. I noticed that when I print zH which is within the function I am finding the root of in both fortran and R I get the value to be in fortran as below

> .\my\_program.exe zH is0.0000000000000000  
> zH is 57.641325192598238  
> zH is 28.820662596299119  
> zH is 14.410331298149559  
> zH is 7.2051656490747797  
> zH is 3.6025828245373899  
> zH is 1.8012914122686949  
> zH is 0.90064570613434747  
> zH is 0.45032285306717373  
> zH is 0.22516142653358687  
> zH is 0.11258071326679343  
> zH is 5.6290356633396717E-002  
> zH is 2.8145178316698358E-002  
> zH is 1.4072589158349179E-002  
> zH is 7.0362945791745896E-003  
> zH is 3.5181472895872948E-003  
> **zH is 1.7590736447936474E-003**  **the last line in fortran**

where as in the R code I got values lower than this

> 1] zH is 0"  
> [1] " zH is 57.6413255382803"  
> [1] " zH is 28.8206627691402"  
> [1] " zH is 14.4103313845701"  
> [1] " zH is 7.20516569228504"  
> [1] " zH is 3.60258284614252"  
> [1] " zH is 1.80129142307126"  
> [1] " zH is 0.90064571153563"  
> [1] " zH is 0.450322855767815"  
> [1] " zH is 0.225161427883908"  
> [1] " zH is 0.112580713941954"  
> [1] " zH is 0.0562903569709769"  
> [1] " zH is 0.0281451784854884"  
> [1] " zH is 0.0140725892427442"  
> [1] " zH is 0.00703629462137211"  
> [1] " zH is 0.00351814731068605"  
> [1] **" zH is 0.00175907365534303"**  **where the fortran ended**  
> [1] " zH is 0.000879536827671514"  
> [1] " zH is 0.000489560168170512"  
> [1] " zH is 0.000674886029990813"  
> [1] " zH is 0.000650298707008735"  
> [1] " zH is 0.00064673969984793"  
> [1] " zH is 0.000646746327972608"  
> [1] " zH is 0.000646746309843868"  
> [1] " zH is 0.000646746309843775"  
> [1] " zH is 0.000646746309843776"  
> [1] " zH is 0.000646746309843776"

After using the one specified in the R code since in R, the default setting is. tol = .Machine$double.eps^0.25.  
I have really been tweaking the ftol value which has been changing the values but stopped at the value above in the call rootscalar\*\*(used 1.0\_wp\*\*(0.25\_wp)\*\* )in fortran. I am working with ecological data and this step is an important one to finish up with the translation.

I would be very grateful if I get suggestions on the next step to take to figure out this as I have been stuck on this step for almost a week. Many thanks in advance

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [September 8, 2023, 3:29pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/2 "2023-09-08T15:29:53Z")

</div>

> [@Mubaraqlanre](#):
>
> … I would be very grateful if I get suggestions on the next step to take to figure out this as I have been stuck on this step for almost a week. Many thanks in advance

@Mubaraqlanre ,

Sorry to read of your trouble, hopefully you will able to resolve it quickly now.

I barely spent a few seconds at your problem and the first glace was at an equation and a thought immediately came to mind you need to apply consistent **precision** between R and Fortran in your arithmetic expressions involving **literal constants** if you want consistent results between the two.

 ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/2X/3/3f6730347c7fe1594f6f8834ec5006a2ec72382b.png)

- Chances are R being effectively an interpreted environment based on GCC ecosystem uses the C language default of **double** for floating-point literal constants.
- Whereas Fortran given its legacy dating back to mid-1950s uses default real which for your compiler might be the equivalent of C **float**.

I suggest you check the precision in R and if it is, as I suspect, the equivalent of C **double** , then reauthor your Fortran code suitably with the use of:

1. defined `KIND` e.g., `use, intrinsic :: iso_c_binding, only : WP => c_double`,

2. use WP as the defined `KIND` consistenly in all floating-point object and constant declarations in Fortran - see point #3 below with `TREF`.

3. use named constants for enhanced mnemonics in your formulae e.g., `real(WP), parameter :: TREF = 298.15_wp` to define the literal constants only once,

4. use the constants such as `TREF` and `INVTREF` in the Fortran equations.

---

<div class="post-metadata">

**Author:** ![Mubaraqlanre](https://avatars.discourse-cdn.com/v4/letter/m/7bcc69/32.png) [@Mubaraqlanre](https://fortran-lang.discourse.group/u/Mubaraqlanre)\
**Post date:** [September 8, 2023, 4:14pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/3 "2023-09-08T16:14:26Z")

</div>

Many thanks for the suggestions and inputs, really appreciate

For 1. I did define the kind and the function is within a subroutine where Tf changing so not a constant. I used the example provided on root fortran sample online. is there need to change this?

> subroutine pH\_calculation  
> use root\_module, only: root\_scalar, wp =\> root\_module\_rk

Regarding the 2. Now I have adjusted the wp instead of kind=8 but I still get the same values

For 3 and 4, I do have parameters defined as the code is part of a bigger model and virtually all values are the same with the R code except for this aspect for the subroutine. For **Tf** in that line was used to refer to a variable within the subroutine that is related to the function, so was not a constant.  
The subroutine is very long but I will share part of it below

> subroutine pH\_calculation  
> use root\_module, only: root\_scalar, wp =\> root\_module\_rk  
> character(len=5), parameter :: method = ‘brent’  
> integer :: ti3, ti  
> real(wp) :: zHV, Tf, V, Vg, x10, x20, x30, y10, y20, y30  
> real (wp) :: y40, InC, InCf, InN, OutN, zHV0, OutC  
> real(wp) :: outnf, InNf, outcf, Zh0f  
> real(wp) :: Tfix, Volumefix, volumegfix, x10f  
> real(wp) :: x20f, x30f, y10f, y20f, y30f, y40f  
> integer:: iflag  
> !real(kind=8), dimension(12) :: roots  
> real(wp) :: result  
> real(wp) :: roots(12)
> 
> ```
> Do ti3=1, cv 
> Do ti = 2,tn
> if (DEBUG_pH) then ! pH module
>     
> if (ti == 2) then
>     
> y10f = 0
>     
> y20f = 0
>     
> y30f = 0
>     
> y40f = 0
>     
> x10f = 0
>     
> x20f = 0
>     
> x30f = 0
> end if
>     
> InCf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0 / 2.0
>                 
> OutNf = (SoilEmissionFlux(ti-1, ti3) * (PatchArea_L/100.0)) / 14.0 + (((TANBudget(ti-1, ti3)) * &
> (1.0-LOSS)) / 14.0)
> 
> InNf = TANproduction(ti, ti3) * (ASL/WSL) / 14.0
> 
> Tfix = Ts(ti) + 273.15 ! K
> 
> Volumefix = WaterBudget2cm(ti, ti3) * (ASL/WSL) ! dm3
> Volumegfix = (Porosity * PatchArea_L * ASL) - Volumefix ! dm3
> 
> OutCf = 0.0
>                
> zH0f = 10.0 ** (-pH(ti-1, ti3)) !Issue with this line
>                 
> 
> Tf=Tfix !To be passed to the function
>                 
> V=Volumefix
>             
> Vg=Volumegfix
>             
> x10=x10f 
>             
> x20=x20f
>             
> x30=x30f
>             
> y10=y10f
>             
> y20=y20f
>             
> y30=y30f
>             
> y40=y40f
>             
> InC=InCf
>             
> InN=InNf
>             
> OutN=OutNf
>             
> zHV0=zH0f
>             
> OutC=OutCf
> !print *, 'Vg is', InCf, ti, ti3
> !print *, 'outN', zH0f, ti, ti3
> call root_scalar('brent',Conc_calc, 0.0_wp, 1.0_wp, zHV, result, iflag, atol=1.0_wp,Ftol = 0.10_wp**(0.25_wp))
>                 
> Hconc(ti, ti3)= zHV ! the root of the function is the H+ concentration
> !print *, 'Vg', ZHV, ti
> pH(ti, ti3)= -log10(Hconc(ti, ti3))
>                  
> roots =Conc_calc2(zHV)
> y10f = roots(4)
> y20f = roots(5)
> y30f = roots(6)
> y40f = roots(7)
> x10f = roots(1)
> x20f = roots(2) 
> x30f =roots(3)
>                 
> !test = reshape([roots(8:12), InC], shape(test))
>               
> else
> pH(ti, ti3) = constpH
> Hconc(ti, ti3) = 10.0 ** (-pH(ti,ti3))
> end if
>        
> end do 
> end do 
> contains 
> function Conc_calc(zHV) result(result)
> real(wp), intent(in) :: zHV
> real(wp) :: A, B, Kw, K1, K2, K3, H1, H2, zH
> real(wp) :: zH0, pH0, pHI, buff_L
> real(wp), dimension(12) :: vec1, Vec2, vec3, vec4, vec5, vec6, vec7
> real(wp), dimension(12) :: vec8, vec9, vec10, vec11, vec12
> real(wp) :: Amatrix(12,12), Bvector(12) 
> real(wp) :: roots(12)
> real(wp) :: result
>         
> integer :: n=12 ! specify the size of the matrix
>     
>     
> A = x10 + x20 + x30 + InN - OutN
>     
> B = y10 + y20 + y30 + y40 + InC
>     
>                
>     
> Kw = 10.0**(-14.0) ! Kw = 10.0**(-14.0) at 25 C
>     
>     
> !Print *, 'zHV', Tf, ti
>                                     
> K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15))) !defines the first H2CO3 dissociation into HCO3-
>     
> K2 = (10.0**(-1.0 * (((2902.39 / Tf) + (0.02379 * Tf) - 6.4980)))) !defines the second HCO3- dissociation into CO3--
>     
> K3 = (10.0**(-1.0 * (((3404.71 / Tf) + (0.032786 * Tf) - 14.8435)))) !defines the third H2CO3-- dissociation into H2O and co2                   
>     
> H1 = (56.0 * exp(4092.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * (1.0 / Tf) * &
>     
> (1.013 * 10.0 ** 5)) 
>     
> H2 = ((1.0 / 29.41) * exp(2400.0 * ((1.0 / Tf) - (1.0 / 298.15)))) / ((1.0 / 8.314) * 0.001 * &
>     
> (1.0 / Tf) * (1.013 * 10.0 ** 5))
> 
>                                       
>     
> zH = zHV * V
>     
> vec1 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, -zH, 0.0_wp, & 
> 0.0_wp, (1.0_wp + H1 * (V / Vg)), 0.0_wp/) ! x1
>         
> if(ti == 2 .and. ti3 == 3) then
> write(*,*) 'vec1 is ', zH             
> end if 
>         
> vec2 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, K1* V, 0.0_wp, 0.0_wp, H1 * (V / Vg), 0.0_wp/) ! x2
> vec6 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/)
> vec3 = (/0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, -zH, 0.0_wp, 0.0_wp, H2* (V / Vg)/) ! y1
>     
> vec4 = (/0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, K2*V, -zH, 0.0_wp, H2 * (V / Vg)/) ! y2
>     
> vec5 = (/1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, &  
> K3* V, 0.0_wp, (1.0_wp + H2 * (V / Vg))/) ! y3
>     
> vec7 = (/0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! y4
>     
> vec8 = (/1.0_wp, 0.0_wp, 0.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH2
>     
> vec9 = (/-1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK3
>     
> vec10 = (/0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK2
>     
> vec11 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RK1
>     
> vec12 = (/0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, -1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp/) ! RH1
>     
>                 
>     
> Amatrix = reshape([vec1, vec2, vec6, vec3, vec4, vec5, vec7, vec8, vec9, vec10, &
>     
> vec11, vec12], shape(Amatrix))
>     
> ! Set the column names
>     
> !character(2) :: colnames(12) = ["x1", "x2", "x3", "y1", "y2", "y3", "y4", "RH2", "RK3", &
>     
> ! "RK2", "RK1", "RH1"]
>     
> Bvector = (/y30, y20 + InC, y10, y40 - OutC, x20 + InN, x10, x30 - OutN, 0.0_wp, 0.0_wp, 0.0_wp, &
> H1 * A * (V / Vg), H2 * B * (V / Vg)/) ! const
>     
> call solve_system(n, Amatrix, Bvector, roots) !I called lapack to solve the system         
>     
> zH0 = zHV0 * V
> pH0 = -log10(zH0 / V)
> pHI = -log10(zH / V) !I renamed the line pH to PHI
> result =(zH0 - InC + (-roots(9) + roots(10) + roots(11)) + (buff_L * (pHI - pH0))) - ((zH))
> end function Conc_calc
> 
> ```
> 
> end subroutine

Thanks a lot for the feedback but could be issue with ftol in root fortran?

---

<div class="post-metadata">

**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\
**Post date:** [September 8, 2023, 5:48pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/4 "2023-09-08T17:48:24Z")

</div>

A recent arXiv preprint [Finding roots of complex analytic functions via generalized colleague matrices](https://arxiv.org/abs/2307.14494) by Zhang and Rokhlin says

> We implemented our algorithms in FORTRAN 77, and compiled it using GNU Fortran,  
> version 12.2.0.

I don’t see the code online, but the authors may be willing to send it if asked.

---

<div class="post-metadata">

**Author:** ![Mubaraqlanre](https://avatars.discourse-cdn.com/v4/letter/m/7bcc69/32.png) [@Mubaraqlanre](https://fortran-lang.discourse.group/u/Mubaraqlanre)\
**Post date:** [September 8, 2023, 7:27pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/5 "2023-09-08T19:27:49Z")

</div>

Thanks a lot for the recommendation.  
I will do a mail to the authors as soon as possible

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [September 8, 2023, 8:22pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/6 "2023-09-08T20:22:27Z")

</div>

> [@Mubaraqlanre](#):
>
> … could be issue with ftol in root fortran?

@Mubaraqlanre ,

I suggest first setting up your formulae **securely** in your Fortran code in a way that will be also **consistent** with the calculations in R.

Consider your equation toward what appears to be enthalpy of second reaction (or is it with Henry’s law constant) that has a lead constant. Your Fortran code has it as `1.0 / 29.41`. You will realize that will result in a difference in the computations in R, likely along the lines shown below:

```fortran
   blk1: block
      real(kind=8) :: a
      a = 1.0 / 29.41
      print *, "Block 1: default precision"
      print "(*(g0))", "a = 1.0/29.41: ", a
      print "(g0,b0)", "a (binary) : ", a
   end block blk1
   print *
   blk2: block
      use, intrinsic :: iso_c_binding, only : WP => c_double
      real(kind=WP), parameter :: ONE = 1.0_wp
      real(kind=WP), parameter :: x = 29.41_wp
      real(kind=WP), parameter :: a = ONE / x
      print *, "Block 2: recommended approach"
      print "(*(g0))", "a = 1.0/29.41: ", a
      print "(g0,b0)", "a (binary) : ", a
   end block blk2
end

```

```auto
C:\temp>gfortran -ffree-form p.f -o p.exe

C:\temp>p.exe
 Block 1: default precision
a = 1.0/29.41: 0.34002039581537247E-1
a (binary) : 11111110100001011010001011011100100000000000000000000000000000

 Block 2: recommended approach
a = 1.0/29.41: 0.34002040122407345E-1
a (binary) : 11111110100001011010001011011100100100101001010110001011001100

```

Note here `block 1` is equivalent your code in Fortran whereas `block 2` is likely equivalent to what occurs in R. Note the difference in the lead constant in the two blocks.

You are dealing with quantities with small numerical magnitudes as pertinent to aqueous ionic chemistry due to limited solubility of CO2 in water leading toward low concentrations of carbonic acid and bicarbonate in the aqueous phase.

So you want to first get the numerical precision consistent in your calculations and then you can proceed toward the solution of the equations and cross-checking the roots you get in Fortran against your R program (the latter albeit at a slower compute speed due to interpreted code at run-time).

And a way to get the numerical precision right in your Fortran code is to get the constants in your equations computed during compilation time with consistently defined `KIND`s as I show in `block 2` above rather than simply writing them as `1.0 / 29.41` in Fortran.

**Caveat emptor with floating-point literal constants in Fortran.**

---

<div class="post-metadata">

**Author:** ![Harper](https://avatars.discourse-cdn.com/v4/letter/h/b5ac83/32.png) [@Harper](https://fortran-lang.discourse.group/u/Harper)\
**Post date:** [September 9, 2023, 12:04am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/7 "2023-09-09T00:04:04Z")

</div>

I saw `Kw = 10.0**(-14.0)` in your program. Why force the machine to do an unnecessary LOG and EXP? If 6-figure accuracy had been enough I would have written `Kw = 1.0e-14` but as you were after higher accuracy I would use `Kw = 1.0e-14_wp`; others have told you about `wp`.

---

<div class="post-metadata">

**Author:** ![gareth](https://avatars.discourse-cdn.com/v4/letter/g/ba8739/32.png) [@gareth](https://fortran-lang.discourse.group/u/gareth)\
**Post date:** [September 9, 2023, 4:23am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/8 "2023-09-09T04:23:45Z")

</div>

As others have stated above, in Fortran real constants are assumed single precision unless you give them a kind suffix – even if the LHS variables are double precision.

So in your code, a line like:

```auto
K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15)))

```

involves lots of single precision calculation on the RHS. So `K1` will have less accuracy than you’d expect from a double precision calculation.

The solution is to provide a suffix every time you use a real constant, e.g.

```auto
K1 = 5.67_wp * (10.0_wp**(-10.0_wp)) * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / 298.15_wp)))

```

In Fortran you should always provide suffixes for real constants, to avoid such problems.

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [September 9, 2023, 7:25am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/9 "2023-09-09T07:25:47Z")

</div>

> [@gareth](#):
>
> The solution is to provide a suffix every time you use a real constant, e.g.
> 
> ```auto
> K1 = 5.67_wp * (10.0_wp**(-10.0_wp)) * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / 298.15_wp)))
> 
> ```

As pointed out elsewhere, the floating point exponents with integer values should almost always be written as integers. However, in this case, there is also the issue of writing the floating point literals in the `usual` way.

```auto
K1 = 5.67e-10_wp * exp(-6286.0_wp * ((1.0_wp / Tf) - (1.0_wp / 298.15_wp)))

```

Of course, it isn’t necessary to do that everywhere because fortran will probably do some of them correctly anyway because of the promotion rules for mixed-kind expressions, but I think it is a good idea for the programmer, especially a new fortran programmer, to specify the kinds explicitly anyway to make his intentions clear.

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [September 10, 2023, 3:16pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/10 "2023-09-10T15:16:29Z")

</div>

Single precision strikes again!

It’s a rite of passage for every Fortran programmer to get bitten by this issue.

TL;DR

- Add `_wp` to every floating point value, no matter how much you don’t want to. The syntax is `3.0_wp`, `10.1e3_wp`, etc.
- Don’t use `3_wp` expecting it to be a floating point value (if it doesn’t have the `.` then it is an integer).
- Don’t use `3.0` without the `_wp` ever under any circumstances.
- Don’t declare variables as `real` or `real(8)` or `real*8` under any circumstances, always `real(wp)`.
- Define `wp`, for example using `use iso_fortran_env, only: wp => real64` to get double precision. Put it in a module if you want to use it in multiple files.

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [September 10, 2023, 11:51pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/11 "2023-09-10T23:51:09Z")

</div>

> [@Mubaraqlanre](#):
>
> … I transitioned a model from R to Fortran due to performance concerns that hindered the efficiency of the original model being to slow in R. … on translating the r code to fortran and solving the code using similar tolerance as used in R code I have not been getting the same values as in the R code. …

@Mubaraqlanre ,

Re: “transitioned a model from R to Fortran due to performance concerns … original model being too slow in R …,” please know you can speed up calculations significantly using Fortran (or, for that matter, C++ or C), basically any compiler based language that processes and optimizes your code directly to your compute platform. Fortran is really good at this once you get over a few of its legacy aspects.

Since your calculation involved the `pH` of the aqueous system of interest containing dissolved CO2, please revisit the system equations (and if you are up to it and can do so, to share them with the readers of this Discourse and any other programming / numerical solution / chemistry forums you follow) and consider:

1. formulating in terms of logarithms (natural log is, well, natural) that extends naturally with `ln K` the natural log of equilibrium constants and the expressions toward them in chemistry literature. You can then solve for the log of hydronium ion concentration(s). This may help you avoid most coding instructions toward exponentiation and the use of `exp` funciton and this can help you further with performance.

2. simple base transformation and/or arithmetic will then yield the `pH` for you.

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [September 11, 2023, 1:04am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/12 "2023-09-11T01:04:40Z")

</div>

So consider the `K` equation you have that “defines the first H2CO3 dissociation into HCO3-” which is in the form of the classic Arrhenius expression:

k = a \ exp(\frac { - E } { RT })  
​  
that is “scaled” in a sense against a reference temperature as

k = a \ exp(\ - E` [\frac { 1 } {T} - \frac { 1 } {T\_{ref}}]\ )  
​  
which you seek to write in Fortran with literal constants as

K1 = 5.67 × 10^{-10} \ exp(\ - 6286.0 \ [\frac { 1 } {T} - \frac { 1 } {298.15}]\ )

So given where things stand with modern Fortran, please see the snippet below that contrasts what you show in the original post with an alternate approach consistent with what R likely calculates for you:

> **Click here for code**
>
> ```fortran
> blk1: block
> real(kind=8) :: K1
> real(kind=8) :: Tf
> print *, "Original post:"
> Tf = 400.0 !
> !defines the first H2CO3 dissociation into HCO3-
> K1 = 5.67 * (10.0**(-10.0)) * exp(-6286.0 * ((1.0 / Tf) - (1.0 / 298.15)))
> print "(*(g0))", "K1(decimal) for Tf = 400 K: ", K1
> print "(g0,b0)", "K1(decimal) for Tf = 400 K: ", K1
> end block blk1
> print *
> blk2: block
> use, intrinsic :: iso_c_binding, only : WP => c_double
> real(kind=WP) :: K1
> real(kind=WP) :: Tf
> real(kind=WP), parameter :: ONE = 1.0_wp
> real(kind=WP), parameter :: TREF = 298.15_wp
> real(kind=WP), parameter :: A1 = 5.67E-10_wp
> real(kind=WP), parameter :: E1 = 6286.0_wp
> print *, "Another approach, more secure in Fortran:"
> Tf = 400.0_wp !
> !defines the first H2CO3 dissociation into HCO3-
> K1 = A1 * exp(-E1 * ( ONE/Tf - ONE/TREF ) )
> print "(*(g0))", "K1(decimal) for Tf = 400 K: ", K1
> print "(g0,b0)", "K1(decimal) for Tf = 400 K: ", K1
> end block blk2
> end
> 
> ```

> **Click here for compiler response**
>
> ```auto
> C:\temp>gfortran -ffree-form p.f -o p.exe
> 
> C:\temp>p.exe
> Original post:
> K1(decimal) for Tf = 400 K: 0.12162608824641047E-6
> K1(decimal) for Tf = 400 K: 11111010000000010100110000101001100010110110100011111110010100
> 
> Another approach, more secure in Fortran:
> K1(decimal) for Tf = 400 K: 0.12162606940969464E-6
> K1(decimal) for Tf = 400 K: 11111010000000010100110000101000111000011011111001111011111000
> 
> ```

Note you are dealing with a weak reaction and the absolute value of the equilibrium constant is small:

- first, notice the difference in the result in the two approaches,
- second, better to work with the natural log of such an equilibrium constant rather than `K1` directly in your equations to be solved.

---

<div class="post-metadata">

**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\
**Post date:** [September 11, 2023, 1:26am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/13 "2023-09-11T01:26:55Z")

</div>

So better yet will be for you to write the scaled expression as

ln\ k = ln\ a\ - (\frac { E` } {T\_{ref}})\ [t\_{r}\ - 1]

where

t\_{r}\ = \frac { T\_{ref} } { T }

which you can author in Fortran as follows and which can help boost your code performance:

> **Click here for code**
>
> ```fortran
> use, intrinsic :: iso_c_binding, only : WP => c_double
> real(kind=WP), parameter :: ONE = 1.0_wp
> real(kind=WP), parameter :: a1 = 5.67E-10_wp ! units?
> real(kind=WP), parameter :: lna1 = log(a1)
> real(kind=WP), parameter :: E1 = 6286.0_wp ! Activation energy, kcal/mol?
> real(kind=WP), parameter :: TREF = 298.15_wp ! degrees Kelvin
> real(kind=WP), parameter :: E1TREF = E1 / TREF
> real(kind=WP) :: K1
> real(kind=WP) :: lnK1
> real(kind=WP) :: Tf
> real(kind=WP) :: Tr
> Tf = 400.0_wp
> Tr = TREF / Tf
> lnK1 = lna1 - E1TREF*( Tr - ONE )
> 
> ```

---

<div class="post-metadata">

**Author:** ![Mubaraqlanre](https://avatars.discourse-cdn.com/v4/letter/m/7bcc69/32.png) [@Mubaraqlanre](https://fortran-lang.discourse.group/u/Mubaraqlanre)\
**Post date:** [September 11, 2023, 7:23am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/14 "2023-09-11T07:23:14Z")

</div>

Thanks to everyone. I have learnt a lot from all the valuable comments and have been adjusting my code appropriately. I have seen some improvements in the value I obtain. I will continue to make changes implementing all the suggestions. Many thanks to everyone for helping with advice.

---

<div class="post-metadata">

**Author:** ![pcosta](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pcosta/32/5504_2.png) [@pcosta](https://fortran-lang.discourse.group/u/pcosta)\
**Post date:** [September 11, 2023, 10:46am UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/15 "2023-09-11T10:46:50Z")

</div>

> [@jacobwilliams](#):
>
> Don’t use `3.0` without the `_wp` ever under any circumstances.

I thought that we won’t face problems when using numbers that can be perfectly described with single precision format, like `3.0`…

```fortran
! cat test.f90
print*,3.d0/acos(-1.d0),3.0/acos(-1.d0)
end

```

Running this program with a recent version of `gfortran` produces the following output:

```shell
0.95492965855137202 0.95492965855137202

```

So in cases like this in mixed precision mode, it is fine, right?

---

<div class="post-metadata">

**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\
**Post date:** [September 11, 2023, 1:13pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/16 "2023-09-11T13:13:45Z")

</div>

Why risk it? Always use `_wp` and you will never have a problem. Also:

- don’t use `.d0` under any circumstances.

---

<div class="post-metadata">

**Author:** ![pcosta](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/pcosta/32/5504_2.png) [@pcosta](https://fortran-lang.discourse.group/u/pcosta)\
**Post date:** [September 11, 2023, 1:23pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/17 "2023-09-11T13:23:14Z")

</div>

Thanks! Though I don’t think I’m risking anything here. I agree with avoiding `.d0`, was just for brevity of the code snippet.

---

<div class="post-metadata">

**Author:** ![ofmla](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ofmla/32/3320_2.png) [@ofmla](https://fortran-lang.discourse.group/u/ofmla)\
**Post date:** [September 11, 2023, 1:23pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/18 "2023-09-11T13:23:19Z")

</div>

within the context of this thread shouldn’t it be `lna1 = dlog(a1)` instead of `lna1 =log(a1)` ? in this case it seems that it doesn’t make difference but I guess that in other cases it does

---

<div class="post-metadata">

**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\
**Post date:** [September 11, 2023, 1:47pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/19 "2023-09-11T13:47:00Z")

</div>

I guess this is true for whole numbers, and in case your processor uses base-2, the set of [dyadic rationals](https://en.wikipedia.org/wiki/Dyadic_rational), i.e. the numbers which can be represented exactly as a fraction of the from i/2^j where i and j are integers. For example `1.5` is 3/2^1 as a fraction.

But on a processor with a different [word-size](https://en.wikipedia.org/wiki/Word_(computer_architecture)), and different floating point format, the constant `1.5` might not be promoted correctly, unless you add the kind specifier suffix.

---

<div class="post-metadata">

**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\
**Post date:** [September 11, 2023, 4:46pm UTC](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489/20 "2023-09-11T16:46:29Z")

</div>

> [@pcosta](#):
>
> I thought that we won’t face problems when using numbers that can be perfectly described with single precision format, like `3.0`…

Yes, and an experienced fortran programmer who understands how the fortran KIND system works, how floating point values of various KINDs are stored, and how mixed-type expressions are converted by default might write code like that. But this discussion is about a new fortran programmer just learning the language, so it is probably better to write out all constants and even all conversions explicitly rather than relying on default conversion conventions.

BTW, there is no difference in the compiled code. Whether you write `x*n` or `x*real(n,kind(x))`, the compiler must do exactly the same thing in both cases, and exactly the same optimizations are allowed. It isn’t like the first expression allows some kind of shortcut that the second expression forbids. And fortran is not an interpreted language, so the number of characters does not affect performance.

[Next page](https://fortran-lang.discourse.group/t/challenge-using-root-fortran-to-find-root-of-a-complex-function/6489.md?page=2)