Please help me understand how to solve the problem "Code of Polar function"

  • Context:
  • Thread starter Thread starter vatlilithuyetvietnam
  • Start date Start date
  • Tags Tags
    Function Polar
Click For Summary

Discussion Overview

The discussion revolves around a coding problem related to the Polar function in Bilayer graphene. Participants are trying to understand why the value of the Polar function becomes excessively large at x=2 and are exploring debugging techniques and potential code modifications to address this issue.

Discussion Character

  • Technical explanation
  • Debugging
  • Mathematical reasoning
  • Exploratory

Main Points Raised

  • One participant reports that the value of the Polar function is large at x=2 and seeks assistance.
  • Another participant suggests stepping through the code to understand its behavior at x=2, implying the importance of debugging.
  • A participant identifies that the denominator of a specific function becomes zero at x=2, indicating a potential source of the problem.
  • One participant attempts to modify the code to handle the case when x=2 to avoid divergence, but expresses frustration that the modification does not work as intended.
  • Another participant recommends wrapping the evaluation of certain functions within conditional statements to avoid calculations that could lead to division by zero.
  • A suggestion is made to avoid exact comparisons with floating-point numbers due to potential precision issues, advocating for a tolerance-based approach instead.
  • One participant reflects on the behavior of the code and acknowledges that certain calculations do not occur when x=2, but also considers the implications of floating-point arithmetic.
  • Another participant emphasizes the importance of documenting any assumptions made regarding the behavior of the code to aid future understanding.

Areas of Agreement / Disagreement

Participants express various viewpoints on how to handle the case when x=2, with some suggesting specific code modifications while others highlight potential pitfalls in floating-point arithmetic. There is no consensus on a definitive solution, and the discussion remains exploratory.

Contextual Notes

Participants note the potential issues with floating-point precision and the need for careful handling of mathematical expressions that could lead to division by zero. There are also concerns about the robustness of the code when dealing with edge cases like x=2.

vatlilithuyetvietnam
Messages
5
Reaction score
0
Hi everybody,
I have a code for Polar function in Bilayer graphene. I compile and run it. Everything is ok but at x=2, the value of Polar function is so large. I don't know why this happens ? Please help me
Code: 
//====================LIBRARY===================================================
#include <stdio.h>
#include <math.h>
#include <conio.h>
//====================CONSTANTS(CGS)============================================
#define  pi               3.141592654
#define  kB               1.3806503*1e-16
#define  hbar             1.0545716*1e-27
#define    msao           0.033*9.10938188*1e-28
#define    elec           4.80320427*1e-10
#define    kappa          2.45
//====================ROUTINE===================================================
double PI(double x,double T1);
//====================MAIN======================================================
double x,T1;
int main()
{
  FILE *output;
  output=fopen ("Pi.dat","w+" );
//   n=2*1e12;
  T1=0.01;
  for (x=0.; x<=4.; x+=0.1)
    {
              printf("%f %f\n ",x,PI(x,T1));
              fprintf(output,"%f %f\n",x,PI(x,T1));
    }
  fclose(output);
  printf("\nfin\n");
  getch();
}
double PI(double x,double T1)
{ double g0,f0,g1,f1,kq;
//  x=2*k1*sin(phi/2);
  g0=sqrt(4+pow(x,4))/2-log((1+sqrt(1+pow(x,4)/4))/2);
  f0=((2+x*x)/2/x)*sqrt(x*x-4)+log((x-sqrt(x*x-4))/(x+sqrt(x*x-4)));
  g1=(1+(pow(x,4)/2)-sqrt(1+pow(x,4)/4))/sqrt(1+pow(x,4)/4);
  f1=((x*x-1)*(pow(x,4)-5*x*x+2))/(x*sqrt(pow(x*x-4,3)));
  if (x>2)
        kq=g0-f0+(pi*pi/6)*(T1*T1)*(g1-f1);
   else if (x<2)
       kq=g0+(pi*pi/6)*(T1*T1)*g1;
   else
          kq= sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-sqrt(2))*(-1.4603545)*sqrt(T1)-sqrt(pi)*(1-sqrt(2)/2)*2.612*pow(sqrt(T1),3);
return kq;
  }
 
Technology news on Phys.org
Well, have you stepped through the code with x=2 to see what it does? Do you understand code debugging?
 
vatlilithuyetvietnam said:
Hi everybody,
I have a code for Polar function in Bilayer graphene. I compile and run it. Everything is ok but at x=2, the value of Polar function is so large. I don't know why this happens ? Please help me

The denominator of
C: 
f1=((x*x-1)*(pow(x,4)-5*x*x+2))/(x*sqrt(pow(x*x-4,3)));
appears to be x(x^2 - 4)^{3/2}. This is zero when x = \pm2. This is a problem.
 
  • Like
Likes   Reactions: Pepper Mint
Thanks phinds, pasmith :)
I tried to test with code:
Code: 
//====================LIBRARY===================================================
#include <stdio.h>
#include <math.h>
#include <conio.h>
//====================CONSTANTS(CGS)============================================
#define  N                1000.
#define  pi               3.141592654
#define  kB               1.3806503*1e-16
#define  hbar             1.0545716*1e-27
#define    msao           0.033*9.10938188*1e-28
#define    elec           4.80320427*1e-10
#define    kappa          2.45
//====================ROUTINE===================================================
double PI(double x,double T1);
double test(double T1);
//====================MAIN======================================================
double x,T1;
int main()
{
  FILE *output;
  output=fopen ("Pi.dat","w+" );
//   n=2*1e12;
  T1=0.01;
  for (x=0.; x<=4.; x+=0.1)
    {
              printf("%f %f %f\n ",x,PI(x,T1),test(T1));
              fprintf(output,"%f %f\n",x,PI(x,T1));
    }
  fclose(output);
  printf("\nfin\n");
  getch();
}
double PI(double x,double T1)
{ double g0,f0,g1,f1,kq;
//  x=2*k1*sin(phi/2);
  g0=sqrt(4+pow(x,4))/2-log((1+sqrt(1+pow(x,4)/4))/2);
  f0=((2+x*x)/2/x)*sqrt(x*x-4)+log((x-sqrt(x*x-4))/(x+sqrt(x*x-4)));
  g1=(1+(pow(x,4)/2)-sqrt(1+pow(x,4)/4))/sqrt(1+pow(x,4)/4);
  f1=((x*x-1)*(pow(x,4)-5*x*x+2))/(x*sqrt(pow(x*x-4,3)));
  if (x==2)
//       kq= sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-sqrt(2))*(-1.4603545)*sqrt(T1)-sqrt(pi)*(1-sqrt(2)/2)*2.612*pow(sqrt(T1),3);
       kq= 1.75486-(0.536077*sqrt(T1))-(1.35599*sqrt(T1)*sqrt(T1)*sqrt(T1));
       else
  {
    if (x<2)
         kq=g0+(pi*pi/6)*(T1*T1)*g1;
   else
       kq=g0-f0+(pi*pi/6)*(T1*T1)*(g1-f1); }
return kq;
  }
double test(double T1)
{ double kq;
  kq= sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-sqrt(2))*(-1.4603545)*sqrt(T1)-sqrt(pi)*(1-sqrt(2)/2)*2.612*pow(sqrt(T1),3);
  return kq;
}

when x=2, I add this line:
Code: 
  if (x==2)
//       kq= sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-sqrt(2))*(-1.4603545)*sqrt(T1)-sqrt(pi)*(1-sqrt(2)/2)*2.612*pow(sqrt(T1),3);
       kq= 1.75486-(0.536077*sqrt(T1))-(1.35599*sqrt(T1)*sqrt(T1)*sqrt(T1));
My purpose is to eliminate divergence at x=2 in my code. But, this command line don't work :( . I don/t know why ? :(
 
f1 is evaluated and used only in case of x!=2, so try wrapping it up in the same scope level of kq. For example
PHP: 
if(x!=2)
{
    if (x>2) //do something, related f,g
    else //do something else , related f,g
}
else //x=2
{
}
Also your loop uses float data x as an iterative value and you make comparison with it in your PI calculation, which is not really good, so you may think about changing your loop method to use integral value as well as your floating point value comparison (we human know 1.9999999!=2.0 but computer might likely not).
Remember to declare and define your objects only within their valid scopes e.g right before they are used, not long before they are used.
 
Last edited:
Thanks you so much, Pepper Mint. Your ideas are so great. I was successful. My code worked.

Code: 
//====================LIBRARY===================================================
#include <stdio.h>
#include <math.h>
#include <conio.h>
//====================CONSTANTS(CGS)============================================
#define  pi               3.141592654
#define  kB               1.3806503*1e-16
#define  hbar             1.0545716*1e-27
#define    msao           0.033*9.10938188*1e-28
#define    elec           4.80320427*1e-10
#define    kappa          2.45
//====================ROUTINE===================================================
double PI(double x,double T1);
//====================MAIN======================================================
double x,T1;
int main()
{
  FILE *output;
  output=fopen ("Pi.dat","w+" );
//   n=2*1e12;
  T1=0.001;
  for (x=0.; x<=4.; x+=0.1)
    {
              printf("%f %f\n ",x,PI(x,T1));
              fprintf(output,"%f %f\n",x,PI(x,T1));
    }
  fclose(output);
  printf("\nfin\n");
  getch();
}
double PI(double x,double T1)
{ double g0,f0,g1,f1,kq;
//  x=2*k1*sin(phi/2);
  g0=sqrt(4+pow(x,4))/2-log((1+sqrt(1+pow(x,4)/4))/2);
  f0=((2+x*x)/2/x)*sqrt(x*x-4)+log((x-sqrt(x*x-4))/(x+sqrt(x*x-4)));
  g1=(1+(pow(x,4)/2)-sqrt(1+pow(x,4)/4))/sqrt(1+pow(x,4)/4);
  f1=((x*x-1)*(pow(x,4)-5*x*x+2))/(x*sqrt(pow(x*x-4,3)));
  if (x>2.0000001)
        kq=g0-f0+(pi*pi/6)*(T1*T1)*(g1-f1);
   else if (x<2)
       kq=g0+(pi*pi/6)*(T1*T1)*g1;
   else
          kq= sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-sqrt(2))*(-1.4603545)*sqrt(T1)-sqrt(pi)*(1-sqrt(2)/2)*2.612*pow(sqrt(T1),3);
return kq;
  }
 
You are doing a lot of calculations of functions of x*x-4 which are very treacherous if the computer thinks that is negative at x=2.0. You should skip all those calculations if x=2.0. But don't test for exact equality. Use a test like (fabs( x-2 ) < 0.000001)
 
Thanks FactChecker due to your ideas :). I think that computer don't see (x*x-4) as a negative number at x=2. Because if x<2, the command line is :
Code: 
else if (x<2)
       kq=g0+(pi*pi/6)*(T1*T1)*g1;
So, f0 doesn't appear, ie sqrt(x*x-4) doesn't appear when we consider case x=2
 
vatlilithuyetvietnam said:
Thanks FactChecker due to your ideas :). I think that computer don't see (x*x-4) as a negative number at x=2. Because if x<2, the command line is :
Code: 
else if (x<2)
       kq=g0+(pi*pi/6)*(T1*T1)*g1;
So, f0 doesn't appear, ie sqrt(x*x-4) doesn't appear when we consider case x=2
You may be right. But unless you want to figure out how the round-off of x*x-4 works inside the processor, I think it is better to be safe. If you figure out that it is safe and execute that code, you should document it so that you and others will not have to figure it out again later. By the time you test x<2 it has already done a lot of powers and sqrts using x*x-4. Are all of those safe?
 
  • #10
I think this has already been alluded to, but f0 and f1 are only used in one if-block, so they should only be calculated in that if-block.
Code: 
   ...
   g0=sqrt(4+pow(x,4))/2-log((1+sqrt(1+pow(x,4)/4))/2);
   g1=(1+(pow(x,4)/2)-sqrt(1+pow(x,4)/4))/sqrt(1+pow(x,4)/4);
   if (x > 2.0)
   {
      double f0=((2+x*x)/2/x)*sqrt(x*x-4)+log((x-sqrt(x*x-4))/(x+sqrt(x*x-4)));
      double f1=((x*x-1)*(pow(x,4)-5*x*x+2))/(x*sqrt(pow(x*x-4,3)));
      kq=g0-f0+(pi*pi/6)*(T1*T1)*(g1-f1);
   }
   else if (x < 2.0)
      kq=g0+(pi*pi/6)*(T1*T1)*g1;
   else
      kq=sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-...
 
  • #11
just a shorter way to merge an ifs...
C: 
if (abs(x-2)<0.000001){
   kq= sqrt(5)-log((1+sqrt(5))/2)-sqrt(pi/4)*(1-...;
}
else {
   double f0=... *(!(x<2.0));
   double f1=... *(!(x<2.0));
   kq= g0-f0+(pi*pi/6)*(T1*T1)*(g1-f1);
}
does that look ok? the if and else if were turned into boolean expressions (that return 0 and 1) since the kq expression looks the same except for the fact that you used the f0,f1. So if x<2.0, !(x<2.0) returns False (or 0) the f0=f1=0 and add nothing to kq (they don't exist) else it's True (or 1) f0=(f0's expression) and f1=(f1's expression) and are included in the kq calculation. They act like Theta's functions.
 

Similar threads

Recursively calculate the standard deviation
  • · Replies 1 ·
Replies
1
Views
2K
Python Finding a local max/min of a function in Python
  • · Replies 7 ·
Replies
7
Views
5K
Problems with Runge-Kutta method
  • · Replies 5 ·
Replies
5
Views
2K
Molecular dynamics Lennard Jones
  • · Replies 4 ·
Replies
4
Views
2K
Verlet algorithm: Why am I getting this output?
  • · Replies 8 ·
Replies
8
Views
4K
C/C++ Trying to access constructor variables in class functions
  • · Replies 36 ·
2
Replies
36
Views
6K
Fortran Verifying Fortran MacCormack Algorithm at t=50*dt
  • · Replies 4 ·
Replies
4
Views
2K
Re-entry Projectile Analysis Program
  • · Replies 7 ·
Replies
7
Views
2K
How to calculate primitive functions on maximal intervals for periodic functions?
  • · Replies 4 ·
Replies
4
Views
2K
Generating Random Numbers with the Acceptance-Rejection Method
  • · Replies 22 ·
Replies
22
Views
5K