Software Projects

Gnu bc truncates intermediate multiply and divide calculations

Posted in Uncategorized by rmt on March 15, 2010

Shell Script to Calculate Protein Radius and Diffusion Distance

#!/bin/bash
#Calculate the radius and diffusion distance of a spherical shape protein
#given the diffusion coefficient, temperature and viscosity of a solution.

pi=3.14159

#Viscosity of water at 298 K and 1 atm (N s m^-2)
v=0.0008931

#Boltzmann constant (J K^-1)
k=1.38065*10^-23

echo "Enter the temperature (K): "
read -e temp
echo "Enter the diffusion coefficient (10^-9 m^2 s^-1): "
read -e d
d=`echo "${d}*10^-09"`
echo "Enter the diffusion time (s): "
read -e dt

r=`echo "scale=30;$k*$temp/($d*6*$pi*$v)" | bc`
dist=`echo "scale=30;sqrt(2*$d*$dt)" | bc`

echo "Radius: $r m" | tee results.txt
echo "Diffusion distance: $dist m" | tee -a results.txt

Despite being an arbitrary-precision language, GNU bc version 1.06 truncates intermediate multiply and divide calculations if the scale= attribute (in lines 21 and 22 here) is not big enough.

Using scale=15 produces an incorrect answer of 0, because 15 is not enough decimal places to hold the Boltzmann constant, which is has a very small coefficient (1.38065×10^-23) that only shows up past the 15th decimal place. Leaving scale= off altogether produces a divide by zero error with the default scale value because the diffusion coefficient gets truncated and treated as a zero in the denominator.

Using scale=30, as is shown above, gives the correct results. (An alternative workaround, to avoid showing all 30 decimal places in the end result, is to compute the numerator and denominator separately using scale=30, and do only the final division step using scale=15.)


Enter the temperature (K):
298
Enter the diffusion coefficient (10^-9 m^2 s^-1):
.069
Enter the diffusion time (s):
60
Radius: .000000003132843654960169968445 m
Diffusion distance: .000090994505328618606651081905 m

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s