Universidad de Panamá
Facultad de ciencias exactas
Alumno:
Bryan Javier Torres Mendoza
Cédula
EC-106-185
Profesor:
Gustavo Jesús Bracho Rodríguez
Materia:
Física computacional
13 de julio de 2022
/* Bryan Javier Torres Mendoza
Física computacional A
Profesor Bracho
Este programa calcula el valor de una integral por el método extendido
del trapecio, el método extendido 3/8 de Simpson y por medio de su
solución exacta, la cual emplea código del libro Numerical Recipes*/
#include <iomanip>
#include <iostream>
#include <cmath>
#include "conio.h"
#include "stdlib.h"
#include <stdio.h>
#include <ctime>
#include <locale.h>
#include "locale.h"
#define ACC 40.0
#include <math.h>
#define EPS 1.0e-10
#define FPMIN 1.0e-30
#define MAXIT 10000
#define XMIN 2.0
#define PI 3.141592653589793
#define NUSE1 5
#define NUSE2 5
#include "nrutil.h"
#define BIGNO 1.0e10
#define BIGNI 1.0e-10
using namespace std;
/* Definición de las variables
b: límite superior.
mu: coeficiente mu.
B: coeficiente beta.
n: número de iteraciones.
x: variable independiente de la función integrada.
*/
///////////////////////////////////////////////////// Función de Chebyshev
float chebev(float a, float b, float c[], int m, float x)
{
void nrerror(char error_text[]);
float d=0.0,dd=0.0,sv,y,y2;
int j;
//if ((x-a)*(x-b) > 0.0) nrerror("x not in range in routine chebev");
y2=2.0*(y=(2.0*x-a-b)/(b-a));
for (j=m-1;j>=1;j--)
{
sv=d;
d=y2*d-dd+c[j];
dd=sv;
}
return y*d-dd+0.5*c[0];
}
/////////////////////////////////////////////////////////// Función beschb
void beschb(double x,double *gam1,double *gam2,double *gampl,double*gammi)
{
float chebev(float a, float b, float c[], int m, float x);
float xx;
static float c1[] = {
-1.142022680371168e0,6.5165112670737e-3,
3.087090173086e-4,-3.4706269649e-6,6.9437664e-9,
3.67795e-11,-1.356e-13};
static float c2[] = {
1.843740587300905e0,-7.68528408447867e-2,
1.2719271366546e-3,-4.9717367042e-6,-3.31261198e-8,
2.423096e-10,-1.702e-13,-1.49e-15};
xx=8.0*x*x-1.0;
*gam1=chebev(-1.0,1.0,c1,NUSE1,xx);
*gam2=chebev(-1.0,1.0,c2,NUSE2,xx);
*gampl= *gam2-x*(*gam1);
*gammi= *gam2+x*(*gam1);
}
/////////////////////////////////////////////////////////// Función bessik
void bessik(float x,float xnu,float *ri,float *rk,float *rip,float *rkp)
{
void beschb(double x, double*gam1, double*gam2,
double*gampl, double*gammi);
void nrerror(char error_text[]);
int i,l,nl;
double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1,gam2,
gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,ripl,
ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2;
//if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessik");
xmu=xnu-nl;
xmu2=xmu*xmu;
xi=1.0/x;
xi2=2.0*xi;
h=xnu*xi;
if (h < FPMIN) h=FPMIN;
b=xi2*xnu;
d=0.0;
c=h;
for (i=1;i<=MAXIT;i++)
{
b += xi2;
d=1.0/(b+d);
c=b+1.0/c;
del=c*d;
h=del*h;
if (fabs(del-1.0) < EPS) break;
}
//if (i > MAXIT) nrerror("x too large in bessik;
//try asymptotic expansion");
ril=FPMIN;
ril1=ril;
rip1=ripl;
fact=xnu*xi;
for (l=nl;l>=1;l--)
{
ritemp=fact*ril+ripl;
fact -= xi;
ripl=fact*ritemp+ril;
ril=ritemp;
}
f=ripl/ril;
if (x < XMIN)
{
x2=0.5*x;
pimu=PI*xmu;
fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu));
d = -log(x2);
e=xmu*d;
fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e);
beschb(xmu,&gam1,&gam2,&gampl,&gammi);
ff=fact*(gam1*cosh(e)+gam2*fact2*d);
sum=ff;
e=exp(e);
p=0.5*e/gampl;
q=0.5/(e*gammi);
c=1.0;
d=x2*x2;
sum1=p;
for (i=1;i<=MAXIT;i++)
{
ff=(i*ff+p+q)/(i*i-xmu2);
c *= (d/i);
p /= (i-xmu);
q /= (i+xmu);
del=c*ff;
sum += del;
del1=c*(p-i*ff);
sum1 += del1;
if (fabs(del) < fabs(sum)*EPS) break;
}
//if (i > MAXIT) nrerror("bessk series failed to converge");
rkmu=sum;
rk1=sum1*xi2;
}
else
{
b=2.0*(1.0+x);
d=1.0/b;
h=delh=d;
q1=0.0;
q2=1.0;
a1=0.25-xmu2;
q=c=a1;
a = -a1;
s=1.0+q*delh;
for (i=2;i<=MAXIT;i++)
{
a -= 2*(i-1);
c = -a*c/i;
qnew=(q1-b*q2)/a;
q1=q2;
q2=qnew;
q += c*qnew;
b += 2.0;
d=1.0/(b+a*d);
delh=(b*d-1.0)*delh;
h += delh;
dels=q*delh;
s += dels;
if (fabs(dels/s) < EPS) break;
}
//if (i > MAXIT) nrerror("bessik: failure to converge in cf2");
h=a1*h;
rkmu=sqrt(PI/(2.0*x))*exp(-x)/s;
rk1=rkmu*(xmu+x+0.5-h)*xi;
}
rkmup=xmu*xi*rkmu-rk1;
rimu=xi/(f*rkmu-rkmup);
*ri=(rimu*ril1)/ril;
*rip=(rimu*rip1)/ril;
for (i=1;i<=nl;i++)
{
rktemp=(xmu+i)*xi2*rk1+rkmu;
rkmu=rk1;
rk1=rktemp;
}
*rk=rkmu;
*rkp=xnu*xi*rkmu-rk1;
}
//////////////////////////////////// Logaritmo natural de la funcion gamma
double gammln(float xx)
{
double x,y,tmp,ser;
static double cof[6]={76.18009172947146,-86.50532032941677,
24.01409824083091,-1.231739572450155,
0.1208650973866179e-2,-0.5395239384953e-5};
int j;
y=x=xx;
tmp=x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.000000000190015;
for (j=0;j<=5;j++) ser += cof[j]/++y;
return -tmp+log(2.5066282746310005*ser/x);
}
////////////////////////////////////////////Solución exacta de la integral
double exacta(double b, double mu, double B)
{
double alfa, x;
float ri, rk, rip, rkp, bessk;
alfa = mu-0.5;
x=(B/(2*b));
bessik(x, alfa, &ri, &rk, &rip, &rkp);
bessk=rk;
return (1/sqrt(PI*b))
* pow(B, 0.5-mu)
* exp(-B/(2*b))
* exp(gammln(mu))
* rk;
}
/////////////////////////////////////////////////////// Función a integrar
double inline funcion(double x, double b, double mu, double B)
{
return pow(x, -2*mu) * pow(b-(x-1e-5), mu-1) * exp(-B/x);
}
/////////////////////////////////////////////////// Derivada de la función
double derivada(double x, double b, double mu, double B)
{
return exp(-B/x)*pow(x,-2*(mu+1))*pow(b-(x-1e-5),mu-2)*(B*(b-x)+x
*(mu*(x-2*b)+x));
}
//////////////////////////////////////////// Error del método del trapecio
double error_trapecio(int n, double a, double b, double mu, double B)
{
return pow(a-b,3) * ((derivada(b,b,mu,B)-derivada(a,b,mu,B))/(b-a))
/ (12*pow(n,2));
}
///////////////////////////////////////////// Regla extendida del trapecio
double Regla_extendida_trapecio(int n,double a,double b,
double mu,double B)
{
double h, suma=0.0;
int i = 1;
h = (b-a)/n;
while (i <= n-1)
{
suma += funcion(a+i*h, b, mu, B);
i++;
}
return (h/2)*(funcion(a, b, mu, B) + 2*suma + funcion(b, b, mu, B))
+ error_trapecio(n, a, b, mu, B);
}
/////////////////////////////////////////// Regla extendida 3/8 de Simpson
double ReglaExt3_8Simpson(double a, double b, int n, double mu, double B)
{
double h, suma1, suma2, suma3, regla38Simpson;
int i, j;
h = (b - a)/(3 * n);
suma1 = suma2 = suma3 = 0.0;
for(i = 1; i <= n - 1; i++) { // Primer sumatorio de la expresión.
j = 3 * i;
suma1 += funcion(a + j * h, b, mu, B);
}
for(i = 1; i <= n; i++) { // Segundo sumatorio de la expresión.
j = 3 * i - 2;
suma2 += funcion(a + j * h, b, mu, B);
}
for(i = 1; i <= n; i++) { // Tercer sumatorio de la expresión.
j = 3 * i - 1;
suma3 += funcion(a + j * h, b, mu, B);
}
regla38Simpson = ((3 * h)/8) * (funcion(a,b,mu,B)+funcion(b,b,mu,B))+
((3 * h)/4) * suma1 + ((9 * h)/8) * suma2 +
((9 * h)/8) * suma3;
//+ Error(a, b, n, h);
return regla38Simpson;
}
//////////////////////////////////////////////////////// Función principal
int main()
{
double b, mu, a, B;
int n;
char opcion;
setlocale(LC_CTYPE,"Spanish");
[Link](6);
[Link](ios::fixed);
system("color F1");
I: cout<<"Introduzca los valores de los coeficientes"<<endl<<endl;
cout<<"Limite superior (b):"; cin>> b;
cout<<"Coeficiente mu (mu):"; cin>> mu;
cout<<"Coeficiente beta (B):"; cin>> B;
cout<<"Cantidad de iteraciones (n): "; cin>> n;
a = 1e-20;
if (a<= 0 || b<=0 || B<=0 || mu<=0)
{
cout<<endl;
cout<<"Error: El valor de a, b y B debe ser un real mayor a 0,"
<<" y el valor de mu debe ser NO menor a 1.";
cout<<endl<<endl<<endl;
goto I;
}
cout<<endl<<"Solución por regla 3/8 extendida de Simpson: "
<<ReglaExt3_8Simpson(a, b, n, mu, B);
cout<<endl<<"Solución por regla extendida del trapecio: "
<<Regla_extendida_trapecio(n, a, b, mu, B);
cout<<endl<<"Solución exacta: "
<<exacta(b, mu, B);
opciones:
{
cout << "\n\n¿Desea volver a realizar cálculos? (S/N): ";
cin >> opcion;
cout << endl;
switch(opcion)
{
case 'S':
case 's':
goto I; // retorna al inicio de la ejecuci?n.
break;
case 'N':
case 'n':
goto final;
break;
default: // Retorna al usuario en caso de opci?n incorrecta.
cout << "\nOpción incorrecta. \n";
cout << "Favor de elegir una de las opciones presentadas. \n";
goto opciones;
}
}
final:
cout<<"\nPrograma finalizado";
return 0;
}
