Integrate - this may be an intermediate question

[agnt1831]

I can't figure out how to NOT get infinity when I integrate through this equation. I need to integrate from 0 to pi/6 but the denominator hits zero at one point.

I have tried a few things to get out of the infinity problem.

Any suggestions? Also, if you have suggestions for my trapezoid rule, I will listen

EDIT: I guess I should have made it more clear. Sorry.

It is 1/(cosx - cosb)


Let's assume b = .523598.... and a = 0 (which is one of the tests I need)

It is 1/(cosx - .866.....)

cosx gets integrated. Cosb remains constant. Since b is one of the bounds of the integral, cosx = cosb at one iteration, which causes the infinity.

I may be missing something more in my coding, but I am not sure.

PS, this is a pendulum integral

//**********TRAPAZOIDAL Pendulum***************




#include <iostream>  
#include <stdio.h>       
#include <math.h> 
using namespace std;


//Modules
//------------------------------------
double f(double x);
double trapezoid(double x);

double b, a, x, dx, xi, t, s, four, xq, tana, t_error, sana, s_error; //I have a few unneded variables here from a previous program
int i, n, q;



int main() {
	

#define PI 3.14159265358979323

	cout.precision(10);

	cout << "Please enter the number of iterations " << endl;
	cin >> n;
	while (n < 0) { //This while loop keeps the iterations greater than zero
		cout << "Please enter a number greater than zero " << endl;
		cin >> n;
	}


	cout << "Please enter an upper bound b " << endl; //Please note, to test PI, enter 3.14159265
	cin >> b;                                         //I defined it above for reference
													  //To keep the code more realistic, I left the bounds as inputs

	cout << "Please enter a lower bound a " << endl;
	cin >> a;

	dx = (b - a) / n;

	cout << "Using the Trapezoidal rule, the integral evaluates to  " << trapezoid(x) << endl;

	




	return dx;


}

//A module with the function that needs to be integrated
//----------------------------------------------
double f(double x) {
	
	

	return 1 / (cos(x) - cos(b));
	

}

//A module to integrate the function using the Trapezoid Rule
double trapezoid(double x) {
	i = 1;
	t = 0;

	while (i > 1 && i < n) {
		xi = (a + (i * dx));
		t = t + (2 * f(xi));
		i = i + 1;
		//cout << xi;
	}

	return ((dx / 2) * (f(b) + t + f(a)));
}

[Ganado]

What function are you trying to integrate using the trapezoidal rule? 1/cos(x)?
Also, there's no reason to have so many global variables. You should define variables in the smallest scope needed, and pass variables to functions as needed.

[jonnin]

I see this: return ((dx / 2) * (f(b) + t + f(a)));

where f(b) is
return 1 / (cos(b) - cos(b)); //because x is b...

I don't know exactly what it should be, but that isnt ever going to work.

[lastchance]

If your function is 1/cos(x) ... which is not what you have coded ... then you should have no difficulty integrating between 0 and pi/6.

What function are you actually trying to integrate?

[agnt1831]

I guess I should have made it more clear. Sorry.

It is 1/(cosx - cosb)


Let's assume b = .523598.... and a = 0 (which is one of the tests I need)

It is 1/(cosx - .866.....)

cosx gets integrated. Cosb remains constant. Since b is one of the bounds of the integral, cosx = cosb at one iteration, which causes the infinity.

I may be missing something more in my coding, but I am not sure.

PS, this is a pendulum integral

[lastchance]

Your pendulum integral is wrong (see * below).

Your integrand should be 1/sqrt(cos x - cos b)

You can avoid the infinity by using the mid-ordinate rule rather than the trapezium rule (or, possibly, transforming the integral).



(*) - pendulum integral:
Kinetic energy + Potential energy = constant
(1/2)mL^2 (d theta / dt)^2 - mgL.cos(theta) = 0 - mgL.cos(b)
which rearranges to give
d theta / dt = sqrt(2g/L)sqrt(cos(theta) - cos(b) )

Then separate variables to get
d theta / sqrt(cos(theta) - cos(b)) = sqrt(2g/L) dt
from which you can integrate.
(x = theta in your code).

[agnt1831]

I can't believe I forgot the sqrt in the integrand. Also, for simplicity, I left the constant out.

Having said that, I need to use the Simpson rule in my code to integrate 1/sqrt(cos x - cos b)

[lastchance]

Personally, I would either:
(a) transform the integral to remove the slightly-badly-behaved integrand; or
(b) use Simpson's rule with either a succession of limits or (slightly more accurately) in conjunction with the mid-ordinate rule at the two ends.

You can't just use Simpson's rule as it stands, it's an improper integral (although it is very-much evaluable).

You could look up "elliptical integrals" if you want to see how to transform your own particular integral.

[agnt1831]

I figured out my problem. There were a few small problems, but the main issue was the last line, where I return. I had an f(b). Obviously, that means the denominator would be ZERO. However, when I removed f(b) from the return, it still gave me bad results. The solution changed with the iterations n.

I switched over to a Simpson Rule calculator I did last week. I didn't want to use this one because I was less certain about my total code.

I removed f(b) and it worked. I replaced f(b) with f(b - .0001). I would get good results at high n if I left f(b) out, but I want it in so I can use it at low n.

It isn't the best looking code, but I think it works. I will go through it once more tomorrow to make sure it's right.

PS, I have a lot of variables I have to remove

//The Pendulum Period Problem using the Simpson Rule

#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;


//Modules
//------------------------------------
double f(double x);
double simpson(double x);

double b, a, x, dx, xi, t, s, four, xq, tana, t_error, sana, s_error;
int i, n, q;



int main() {

#define PI 3.14159265358979323
#define g 9.81
#define l 9.81

cout.precision(10);

cout << "Please enter the number of iterations " << endl;
cin >> n;
while (n < 0) { //This while loop keeps the iterations greater than zero
cout << "Please enter a number greater than zero " << endl;
cin >> n;
}

cout << "Please enter an upper bound b " << endl;
cin >> b;

cout << "Please enter a lower bound a " << endl;
cin >> a;

dx = (b - a) / n;


cout << "The period of the pendulum is " << sqrt(8*l/g) * (simpson(x) + ((dx / 3) * f(b-.0001))) << endl;

return dx;
}

//A module to calculate the original function
//----------------------------------------------
double f(double x) {

//*****Please insert new equation******
return (1 / sqrt(cos(x) - cos(b))); // <--------------
}




double simpson(double x) {
q = 1; // counter
s = 0; //equation
four = 4;


while (q >= 0 && q < (n-1)) {
xq = a + (q * dx);
s = s + ((dx / 3) * four * f(xq));
q = q + 1;
if (q % 2 == 0) {

four = 2;
}
else {
four = 4;

}
}

return s + ((dx / 3) * f(a));
}

[lastchance]

You would do far better to transform the integral to avoid the singularities.

Use substitution

cos(x)-cos(b) = 2 ( sin(b/2) cos(u) )^2

The integrand transforms to

sqrt( 2 / ( 1 - (sin(b/2)sin(u))^2 ) )

and the limits to 0 and pi/2. There is no longer a singularity in the integral.


For small angles you can compare it with the simple-harmonic-motion (SHM) approximation.

#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;

const double PI = 4.0 * atan( 1.0 );
double beta;                           // no easy way round this global variable


double Simpson( double a, double b, double (*f)( double ), int n )
{
   if ( n < 4 || n % 2 ) { cerr << "Invalid n";    exit(1); }

   double dx = ( b - a ) / n;
   double sumOdd = 0.0, sumEven = 0.0;
   for ( int i = 1; i < n    ; i +=2 ) sumOdd  += f( a + i * dx );
   for ( int i = 2; i < n - 1; i +=2 ) sumEven += f( a + i * dx );
   return ( f( a ) + f( b ) + 4 * sumOdd + 2 * sumEven ) * dx / 3;
}


double fn( double theta )
{
   double s = sin( beta / 2.0 ) * sin( theta );
   return 1.0 / sqrt( 1.0 - s * s );
}


int main()
{
   const double g = 9.81;   // gravity (m/s^2)
   const double L = 1.0;    // pendulum length (m)

   int n;
   cout << "Enter number of intervals (even, greater than 2): ";   cin >> n;
   cout << "Enter maximum angular displacement (radians): ";   cin >> beta;

   cout << "Period of pendulum = " << 4.0 * sqrt( L / g ) * Simpson( 0.0, 0.5 * PI, fn, n ) << '\n';
   cout << "Small-angle approximation (SHM): " << 2.0 * PI * sqrt( L / g ) << '\n';
}

Small-ish angle displacement:

Enter number of intervals (even, greater than 2): 100
Enter maximum angular displacement (radians): 0.1
Period of pendulum = 2.00732
Small-angle approximation (SHM): 2.00607

pi/6 angle displacement

Enter number of intervals (even, greater than 2): 100
Enter maximum angular displacement (radians): 0.5236
Period of pendulum = 2.04099
Small-angle approximation (SHM): 2.00607