Implementation of pow function in c language with power of -3/2 or -5/2

Question

I am new to C language. I am trying to implement the formula (r to power of (-3/2)) * exp(t*(r-1/r). I am using pow function but when I keep this equation in for loop it continuously implementing the loop saying that "pow domain error". I want the increment in for loop to be of 0.2. Please help me in implementing this formula. The code is:

#include <stdio.h>
#include <conio.h>
#include <math.h>

void main()
{
    double R,r;
    float f;
    clrscr();
    for(r=0;r<=5;r+0.2)
    {
        R=pow(r,-1.5)*exp(2.5*(r-1)/r);
        f=r-R;
        printf("value of f is:%f",f);
    }
    getch();
}

Change `for(r=0;r<=5;r+0.2)` to `for(r=1;r<=5;r+=0.2)`. (The computation of `pow(r,-1.5)` and `(r-1)/r` will fail when r = 0.) — r3mainer, Mar 15 '18 at 22:12
Your loop `for(r=0;r<=5;r+=0.2)` will end with iteration in which `r` is approximately 4.8; it will not execute an iteration `r` is 5 (or approximately 5). This is in part because when .2 is converted to the floating-point format common C implementations use, the result is 0.200000000000000011102230246251565404236316680908203125. To iterate properly, use `for (double i = 0; i <= 50; i += 2)` and put `r = i / 10;` inside the loop. — Eric Postpischil, Mar 15 '18 at 23:04

sg7 · Accepted Answer · 2018-03-16T12:16:19.487

You may want to be aware that:

The C Standard, 7.12.1 [ISO/IEC 9899:2011], defines three types of errors that relate specifically to math functions in . Paragraph 2 states

A domain error occurs if an input argument is outside the domain over which the mathematical function is defined.

Paragraph 3 states

A pole error (also known as a singularity or infinitary) occurs if the mathematical function has an exact infinite result as the finite input argument(s) are approached in the limit.

Paragraph 4 states

A range error occurs if the mathematical result of the function cannot be represented in an object of the specified type, due to extreme magnitude.

The available domain for function pow(x,y) is:

 (x > 0) || (x == 0 && y > 0) || (x < 0 && y is an integer)

Which means that you cannot use r = 0 as an argument of pow.

1) pow(r, -1.5) = r^(-3/2) = 1.0/ r^3/2 = 1.0/sqrt(r*r*r);

2) In C language main function should return int value and takes void not empty bracket ()

3) Incrementing operation can be done via r = r + 0.2 or r += 0.2

4) To avoid division by zero in (r-1)/r close to 0 value has to be used instead.

5) In order to not loose double precision the variable f should be also declared double.

6) The value of f for r=0 is extremely difficult to calculate. The algorithm is numerically unstable for r = 0. We have division of two infinities. The proper way to solve it is to find function limit at point 0+. Numerically, its is enough to set r close to 0.

7) You may consider much more smaller step for r in order to capture a better plot due to the nature of your f function.

#include<stdio.h>
#include<conio.h>
#include<math.h>

#define  ZERO_PLUS      1e-14       

void calculate_and_print(double r)
{
    double f, f1;

    f =  r -       pow(r, -1.5) * exp(2.5*(r-1)/r);
    f1 = r - (1.0 /sqrt(r*r*r)) * exp(2.5*(r-1)/r);

    printf("Value of f for r= %f  is  f= %f  f1= %f \n", r, f, f1);
}

int main(void)
{
    clrscr();

    calculate_and_print(ZERO_PLUS);  

    for (double r = 0.2; r < 5.0; r += 0.2 )
        calculate_and_print(r);

    calculate_and_print(5.0); 

    return 0;
}

Test:

Value of f for r= 0.000000  is  f= 0.000000  f1= 0.000000                                                                                     
Value of f for r= 0.200000  is  f= 0.199492  f1= 0.199492                                                                                     
Value of f for r= 0.400000  is  f= 0.307038  f1= 0.307038                                                                                     
Value of f for r= 0.600000  is  f= 0.193604  f1= 0.193604                                                                                     
Value of f for r= 0.800000  is  f= 0.051949  f1= 0.051949                                                                                     
Value of f for r= 1.000000  is  f= 0.000000  f1= 0.000000                                                                                     
Value of f for r= 1.200000  is  f= 0.046058  f1= 0.046058                                                                                     
Value of f for r= 1.400000  is  f= 0.166843  f1= 0.166843                                                                                     
Value of f for r= 1.600000  is  f= 0.338256  f1= 0.338256                                                                                     
Value of f for r= 1.800000  is  f= 0.542116  f1= 0.542116                                                                                     
Value of f for r= 2.000000  is  f= 0.765977  f1= 0.765977                                                                                     
Value of f for r= 2.200000  is  f= 1.001644  f1= 1.001644                                                                                     
Value of f for r= 2.400000  is  f= 1.243810  f1= 1.243810                                                                                     
Value of f for r= 2.600000  is  f= 1.489073  f1= 1.489073                                                                                     
Value of f for r= 2.800000  is  f= 1.735278  f1= 1.735278                                                                                     
Value of f for r= 3.000000  is  f= 1.981075  f1= 1.981075                                                                                     
Value of f for r= 3.200000  is  f= 2.225642  f1= 2.225642                                                                                     
Value of f for r= 3.400000  is  f= 2.468498  f1= 2.468498                                                                                     
Value of f for r= 3.600000  is  f= 2.709387  f1= 2.709387                                                                                     
Value of f for r= 3.800000  is  f= 2.948194  f1= 2.948194                                                                                     
Value of f for r= 4.000000  is  f= 3.184898  f1= 3.184898                                                                                     
Value of f for r= 4.200000  is  f= 3.419535  f1= 3.419535                                                                                     
Value of f for r= 4.400000  is  f= 3.652177  f1= 3.652177                                                                                     
Value of f for r= 4.600000  is  f= 3.882916  f1= 3.882916                                                                                     
Value of f for r= 4.800000  is  f= 4.111856  f1= 4.111856                                                                                     
Value of f for r= 5.000000  is  f= 4.339103  f1= 4.339103

Let me know if you have more questions.

One can loop with integer values in `double` format. That is, they ought to increment by 2 instead of .2 but may do it in a floating-point type instead of an integer type. Floating-point handles integers exactly, up to the point where they overflow the width available (same as integer formats), and using floating-point instead of integer avoids repeated integer-to-floating-point conversions. — Eric Postpischil, Mar 16 '18 at 01:39
@EricPostpischil You are absolutely right! I corrected the answer and increment the double to avoid conversions. I increment by `1.0` in order to avoid division by `10.0`. The multiplication is faster then the division. — sg7, Mar 16 '18 at 01:56
how can we plot the graph between f and r using c language .Is there any function that helps us plotting . — Sudhir Kumar, Mar 16 '18 at 19:53
@SudhirKumar C language does not have any build in plot function. An external library would have to be used. Check: https://stackoverflow.com/questions/1275484/good-plotting-library-for-c — sg7, Mar 16 '18 at 22:41
when i am trying to write these values to a file it says the pointer is not defined and just only one values is getting printed to file — Sudhir Kumar, Mar 22 '18 at 15:36
@SudhirKumar I will be happy to review your code. Post it as a separate questions. Do you need more help with the current question? — sg7, Apr 01 '18 at 01:28

score 1 · Answer 2 · answered Mar 15 '18 at 22:22

1

From http://en.cppreference.com/w/c/numeric/math/pow:

If base is zero and exp is negative, a domain error or a pole error may occur.

You run into that situation when you evaulate pow(r, -1.5) with the value of r set to zero.

Change your loop so that r starts with 0.2 instead of zero.

Also, you need to use r += 0.2 instead of r+0.2 in the for statement.

for ( r = 0.2; r <= 5; r += 0.2)
{
   R = pow(r, -1.5)*exp(2.5*(r-1)/r);
   f = r-R;
   printf("value of f is:%f", f);
}

answered Mar 15 '18 at 22:22

R Sahu

204,454
14
159
270

once i get the results i should be plotting the graph for r and f.So the initial value that r has to be 0 . – Sudhir Kumar Mar 15 '18 at 23:12
@SudhirKumar Well, if you can figure out what `pow(0, -1.5)` is, you can then move on to the problem of what `2.5*(0-1)/0` is. – aschepler Mar 16 '18 at 01:48

Implementation of pow function in c language with power of -3/2 or -5/2

2 Answers2