Optimizing Electric Potential Calculation for a Charged Particle in 1-Dim

[Saladsamurai]

so the following program uses a user defined function called epotential which is basically a function of the form f(x)=c/x where c is a constant.

I am using a for loop to evaluate the function from -0.5 to +0.5 and incrementing by 0.1

The problem is that I do not want it to evalute f(0) since it is undefined.

I thought that I could simply put "if (x!=0)...evaluate the function"

but this is not working and I do not see what is wrong with that logic.

Can someone help me out here?

#include <iostream>
#include <fstream>         // needed for outfile
using namespace std;


double epotential (double charge, double dist);     // function prototype

int main ()

{
    ofstream outfile;             
    outfile.open ("ePotential.txt");           // opens textfile
    
    //declare and initialize objects
    
    double charge=8.0*0.000000001;
    double dist;
    double V;
    
    
    outfile <<"x= " << '\t' << "V= " <<endl << endl;         //cout to textfile
    cout << "x= " << '\t' << "V= " <<endl << endl;           //cout to screen
    
    double inc = 0.1;            // increment in for loop
    double left_bound = -0.5;    // left bound of domain of x
    double right_bound = 0.5;    // right bound of domain of x
    

    for (dist= left_bound; dist <= right_bound; dist+=inc)
    
    {
           if (dist != 0)  // I do not want it to output V if dist = 0
                           // why does this not work?!
                           
           V=epotential(charge, dist);               // loops over domain of x
           outfile << dist << '\t' << V << endl;     // sends output to textfile   
           cout << dist << '\t' << V << endl;        // sends output to screen
        
    }
    
    system ("PAUSE");                    // for Dev C++ compiler
    return 0;
}
        
        


        
    
    
//function that calculates the electric potential of a charged particle in 1-dim
/*---------------------------------------------------------------------*/
/*---------------------------------------------------------------------*/

double epotential (double q, double r)  
{
       double V;
       double k=9.0*1000000000;
       
       V=k*(q/r);
       
       return (V);
}

Here is the output

x= V=

-0.5 -144
-0.4 -180
-0.3 -240
-0.2 -360
-0.1 -720
-2.77556e-017 -2.59407e+018
0.1 720
0.2 360
0.3 240
0.4 180
0.5 144

 

[rootX]

if (dist != 0) // I do not want it to output V if dist = 0
// why does this not work?!

V=epotential(charge, dist); // loops over domain of x
outfile << dist << '\t' << V << endl; // sends output to textfile
count << dist << '\t' << V << endl; // sends output to screen

How the compiler knows when to stop if?

Or, how would you write the following code?

int main ()
{
if condition met:
do something

do something all times
}

 

[mgb_phys]

You can't compare the rest of a floating point calculation to 0
As dist is incremented it might equal 0.00000000000001 an so the test fails
You need to convert to an int or check for a range.

 

[Saladsamurai]

How the compiler knows when to stop if?

Or, how would you write the following code?

int main ()
{
if condition met:
do something

do something all times
}

I am sorry I don't follow. I more or less have the following:

for x=-0.5 to 0.5 (increment by 0.1)

if x not equal to 0
evaluate V(x)

You can't compare the rest of a floating point calculation to 0
As dist is incremented it might equal 0.00000000000001 an so the test fails
You need to convert to an int or check for a range.

Sorry don't follow this either? I am an amateur, so bare with me. I am looping through the domain -0.5 to 0.5 and incrementing by 0.1 ... so at one point dist = 0 .

So what is the problem here?

 

[rootX]

I am sorry I don't follow. I more or less have the following:

for x=-0.5 to 0.5 (increment by 0.1)

if x not equal to 0
evaluate V(x)

Forgetting the brackets.

V=epotential(charge, dist); // loops over domain of
outfile << dist << '\t' << V << endl; // sends output to
count << dist << '\t' << V << endl; // sen

In my VS, if I don't include brackets, only first line below the if is considered a part of it. You have three lines and you define V to be under if, but then accessing V outside the scope of if condition. I only read that part of the code but wonder if this would create any problems

Sorry don't follow this either? I am an amateur, so bare with me. I am looping through the domain -0.5 to 0.5 and incrementing by 0.1 ... so at one point dist = 0

There's nothing like exact 0 , 0.5, or 0.1.

 

[Saladsamurai]

Forgetting the brackets.

V=epotential(charge, dist); // loops over domain of
outfile << dist << '\t' << V << endl; // sends output to
count << dist << '\t' << V << endl; // sen

In my VS, if I don't include brackets, only first line below the if is considered a part of it. You have three lines and you define V to be under if, but then accessing V outside the scope of if condition. I only read that part of the code but wonder if this would create any problems

Okay, I think that I tried that. I understand now.

There's nothing like exact 0 , 0.5, or 0.1.

Do you mean when it adds 0.1 + -0.1 it does not equal EXACTLY zero ?

Is that what you are saying?

 

[rootX]

Do you mean when it adds 0.1 + -0.1 it does not equal EXACTLY zero ?

Is that what you are saying?

Yes. See ~third post:

As dist is incremented it might equal 0.00000000000001 an so the test fails

http://idlastro.gsfc.nasa.gov/idl_html_help/Accuracy_and_Floating_Point_Operations.html

 

[Saladsamurai]

Yes. See ~third post:


http://idlastro.gsfc.nasa.gov/idl_html_help/Accuracy_and_Floating_Point_Operations.html

Gotcha. Thanks!

I am not sure, though, how I can convert this to int ...?

 

[mgb_phys]

Either do
if ((int)dist == 0 ) or
if ( (dist > -0.00001) && (dist < 0.00001) )

 

[Saladsamurai]

Either do
if ((int)dist == 0 ) or
if ( (dist > -0.00001) && (dist < 0.00001) )

You mean if((int)dist != 0) right?

Since I do not want f(0)

 

[Saladsamurai]

This is getting silly.

If I replace if (dist != 0)with if ((int)dist != 0) nothing from the loop gets outputted
if I replace it with

if ((int)dist == 0) i get the same output as before

Either do
if ((int)dist == 0 ) or
if ( (dist > -0.00001) && (dist < 0.00001) )

if I want to use either of these, I have to tell the compiler that if (int)dist == 0

then DONT evaluate f(dist)

and I cannot figure out how to do that. Any ideas?

 

[Saladsamurai]

Ha ha! Gotcha!

 for (dist= left_bound; dist <= right_bound; dist+=inc)
    
    {
           if ((dist > -0.00001) && (dist < 0.00001))
           {continue;}
                           
           else                
            {
            V=epotential(charge, dist);               // loops over domain of x
            outfile << dist << '\t' << V << endl;     // sends output to textfile   
            cout << dist << '\t' << V << endl;        // sends output to screen
           }      
    }

There has to be a way using int though...

 

[rootX]

Ha ha! Gotcha!

 for (dist= left_bound; dist <= right_bound; dist+=inc)
    
    {
           if ((dist > -0.00001) && (dist < 0.00001))
           {continue;}
                           
           else                
            {
            V=epotential(charge, dist);               // loops over domain of x
            outfile << dist << '\t' << V << endl;     // sends output to textfile   
            cout << dist << '\t' << V << endl;        // sends output to screen
           }      
    }

There has to be a way using int though...

Well, int converts float to an integer: -1,0,1
so int (0.5) would obviously round down to 0 not 1.
I think one of the way to use int would be multiply by like 1000 before converting to int.

 

[D H]

You don't need to use integers.

As you saw, two numbers that are the same on paper can be different on a computer. It is usually not a good idea to compare two doubles for equality. It is much safer to test for "reasonably close to the same", where "reasonably close" application specific. Here you are incrementing by 0.1 (more on this below) from -0.5, so the fifth increment should put you somewhere close to zero. So test for that -- zero plus or minus some small epsilon.

The problem here is that 0.1 (decimal) in binary is 0.000110001100011... It's binary representation is non-terminating, just as 1/3 has a non-terminating decimal expansion. A rational number with a non-terminating binary expansion cannot be represented exactly using the standard computer representation for doubles. By comparing with ==, you are acting as if the computation is exact. It isn't.

 

[Coin]

Well, int converts float to an integer: -1,0,1
so int (0.5) would obviously round down to 0 not 1.
I think one of the way to use int would be multiply by like 1000 before converting to int.

This would also be a good time to pull out the floor(), ceiling() and round() functions from math.h.

(Although come to think of it is it guaranteed that round( n ) == 0 where -0.5 < n < 0.5 ?)

 

[AUMathTutor]

You could always try a trick like this:

for(int i = 0; i < 11; i++)
{
   if(i != 0)
   {
      double dist = (i - 5) * 0.1;
      // code...
   }
}

Essentially, since you will be using a finite discrete number of values for dist, there exists a mapping from integers to reals. And since integers are computed exactly (in the ranges, anyway) this will work *every* time. You could even generalize the code to work with arbitrary bounds. For instance...

double proportionalMapping(double min, double max, int sampleSize, int i)
{
   return min + (max - min) * i / (sampleSize-1);
}

Then your original code could look like

const int sampleSize = 11;
const double min = -0.5, max = 0.5;

for(int i = 0; i < sampleSize; i++)
{
   if(i != 0)
   {
      double dist = proportionalMapping(min, max, sampleSize, i);
   }
}