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

In summary, the conversation discusses a code for a Polar function in Bilayer graphene and the issue of a large value appearing at x=2. The code is compiled and run, but the problem persists. The conversation also includes suggestions for debugging the code and modifying the loop method to avoid floating point comparisons. The final code is shown to have successfully addressed the problem.
  • #1
vatlilithuyetvietnam
5
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
  • #2
Well, have you stepped through the code with x=2 to see what it does? Do you understand code debugging?
 
  • #3
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 [itex]x(x^2 - 4)^{3/2}[/itex]. This is zero when [itex]x = \pm2[/itex]. This is a problem.
 
  • Like
Likes Pepper Mint
  • #4
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 ? :(
 
  • #5
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:
  • #6
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;
  }
 
  • #7
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)
 
  • #8
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
 
  • #9
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.
 

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

1. What is a Code of Polar function?

A Code of Polar function is a mathematical representation of a polar coordinate system, where a point on a plane is defined by its distance from the origin and its angle from a fixed reference line. In coding, this function is often used to convert between polar and Cartesian coordinates.

2. How do I solve a Code of Polar function?

To solve a Code of Polar function, you will need to use mathematical formulas and equations. First, you will need to understand the basics of polar coordinates and how they relate to Cartesian coordinates. Then, you can use trigonometric functions and the Pythagorean theorem to convert between the two coordinate systems and solve the problem.

3. What are some common applications of Code of Polar function?

Code of Polar function is commonly used in computer graphics, navigation systems, and robotics. It is also used in physics and engineering to describe circular or rotational motion.

4. Can you provide an example of how to solve a Code of Polar function?

Sure, let's say we have a polar function r = 3cos(theta). To convert this into Cartesian coordinates, we can use the equations x = rcos(theta) and y = rsin(theta). So, plugging in 3cos(theta) for r, we get x = 3cos^2(theta) and y = 3cos(theta)sin(theta). This gives us a graph of a circle with a radius of 3 centered at the origin.

5. Are there any tips or tricks for solving Code of Polar function problems?

One helpful tip is to always start by drawing a diagram to visualize the problem. This can make it easier to understand the relationship between the polar and Cartesian coordinates. Also, make sure to review the formulas and equations used to convert between the two coordinate systems. Practice and repetition can also improve your understanding and problem-solving skills for Code of Polar function.

Similar threads

Monte-Carlo integration does not work properly after implementing pdfs
Replies
1
Views
816
Recursively calculate the standard deviation
Replies
1
Views
2K
Problems with Runge-Kutta method
Replies
5
Views
2K
Python Finding a local max/min of a function in Python
Replies
7
Views
4K
Molecular dynamics Lennard Jones
Replies
4
Views
1K
Verlet algorithm: Why am I getting this output?
Replies
8
Views
3K
Fortran Verifying Fortran MacCormack Algorithm at t=50*dt
Replies
4
Views
2K
Re-entry Projectile Analysis Program
Replies
7
Views
1K
Trying to compile a source code in C
Replies
5
Views
2K
MATLAB MATLAB IFFT doesn't match the analytical one
Replies
10
Views
2K

Hot Threads

Recent Insights

Back
Top