I don't know how to implement a function.
I have a lot of trouble when dealing with this problem.
I would be rather happy if anyone could explain it to me.
/*
##############################################################################
FUNCIÓN PARA CALCULAR LA DERIVADA PRIMERA DE UNA FUNCIÓN USANDO LAS FÓRMULAS VISTAS EN CLASE
- LOS PUNTOS VIENEN DADOS EN EL VECTOR x.
- LA FUNCIÓN f(x) VIENE DADA POR EL VECTOR fx,
- HAY QUE TENER EN CUENTA LA GESTIÓN DE LOS BORDES, ES DECIR, DE LOS VALORES DE LOS EXTREMOS DEL VECTOR:
- BORDES: CALCULAR DERIVADA CON 2 PUNTOS. ES DECIR : f'(xi)=(fx(xj)-fx(xi))/(xj-xi)
- RESTO DE ELEMENTOS: CALCULAR DERIVADA CON 3 PUNTOS. ES DECIR:
f'(xi)=((xi-xl)(f(xr)-f(xi))/(xr-xi) + (xr-xi)(f(xi)-f(xl))/(xi-xl)) / (xr-xl)
- LA DERIVADA SE DEVUELVE EN EL VECTOR DE ENTRADA/SALIDA fd.
- LA FUNCIÓN DEVUELVE 0 SI TERMINA CORRECTAMENTE, Y -1 EN CASO DE ERROR.
###############################################################################
*/
// #############################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR EL MÉTODO DEL TRAPECIO
// UTILIZANDO N INTERVALOS EN UN INTERVALO [a,b]
// #############################################################
real mn_trapecio(
real (*f)(const real x), // Función que se integra
const real a, // Extremo izquierdo del intervalo
const real b, // Extremo derecho del intervalo
const int N) // Número de intervalos para calcular la integral
{
//HACER ALUMNO
real integral = 0; // Variable donde se acumula el valor de la integral
real funcion = (*f)(a); // Evaluación de la función en a
real h = (b-a)/(2.*N); // Factor multiplicativo
for (int i = 0; i < N; i++) { // Bucle para acumular el valor de la integral
real ext_dch = a + (i+1) * (b-a)/N; // Extremo derecho subintervalo
real ev_funcion = (*f)(ext_dch); // Evaluación de la función en el extremo derecho
integral += h * (funcion + ev_funcion); // Acumulación valor integral
funcion = ev_funcion; // Actualización valor integral extremo izquierdo
}
return integral; // Se devuelve el valor de la integral
}
// ###########################################################################
// FUNCION PARA CALCULAR UNA INTEGRAL POR ITERACIONES DEL METODO DEL TRAPECIO
// SE PARTE DE N=1 INTERVALOS Y SE CALCULA LA INTEGRAL CON LA FUNCION ANTERIOR
// SE ACTUALIZA N MULTIPLICANDOLO POR 2 Y SE CALCULA DE NUEVO LA INTEGRAL
// SE REALIZA ESTE PROCESO HASTA QUE LA DISTANCIA RELATIVA ENTRE LOS 2
// VALORES DE LA INTEGRAL ES MENOR QUE LA TOLERANCIA. EN Niter SE DEVUELVE
// EL NUMERO DE ITERACIONES QUE SE REALIZA DE ESTE PROCESO
// ################################################################
mn_trapecio(
real (*f)(const real x), // Función que se integra
const real a, // Extremo izquierda del intervalo
const real b, // Extremo derecho del intervalo
int &Niter, // Variable de salida con el número de iteraciones realizadas
const real TOL) // Tolerancia para controlar la convergencia
{
// HACER ALUMNO
int N = 1;
real integral = mn_trapecio((*f),a,b,N);
real error = TOL + 1;
Niter = 0;
while (error > TOL) {
if(N > 1e7) break; // Fijamos "a priori" el número máximo de intervalos a 1e7
Niter++;
N *= 2;
real integral2 = mn_trapecio((*f),a,b,N);
error = mn_distancia(integral, integral2);
integral = integral2;
}
return (integral);
}
// ###########################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR EL MÉTODO DE SIMPSON
// ###########################################################
real mn_simpson(
real (*f)(const real x), // Función que se integral
const real a, // Extremo izquierdo del intervalo
const real b, // Extremo derecho del intervalo
const int N) // Número de intervalos para calcular la integral
{
//HACER ALUMNO
real integral = 0; // Variable donde se acumula el valor de la integral
real funcion = (*f)(a); // Evaluación de la función en a
real h = (b-a)/(6.*N); // Factor multiplicativo
for (int i = 0; i < N; i++) { // Bucle para acumular el valor de la integral
real ext_izq = a + i * (b-a)/N; // Extremo izquierdo subintervalo
real ext_dch = a + (i+1) * (b-a)/N; // Extremo derecho subintervalo
real ev_funcion = (*f)(ext_dch); // Evaluación de la función en el extremo derecho
integral += h * (funcion + ev_funcion + 4 * (*f)((ext_izq + ext_dch)/2.)); // Acumulación valor integral
funcion = ev_funcion; // Actualización valor integral extremo izquierdo
}
return integral; // Se devuelve el valor de la integral
}
// ########################################################################
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR ITERACIONES DEL MÉTODO DE SIMPSON
// ########################################################################
real mn_simpson(
real (*f)(const real x), // Función que se integra
const real a, // Extremo izquierda del intervalo
const real b, // Extremo derecho del intervalo
int &Niter, // Variable de salida con el número de iteraciones realizadas
const real TOL) // Tolerancia para controlar la convergencia
{
// HACER ALUMNO
int N = 1;
real integral = mn_simpson((*f),a,b,N);
real error = TOL + 1;
Niter = 0;
while (error > TOL) {
if(N > 1e7) break; // Fijamos "a priori" el número máximo de intervalos a 1e7
Niter++;
N *= 2;
real integral2 = mn_simpson((*f),a,b,N);
error = mn_distancia(integral, integral2);
integral = integral2;
}
return (integral);
}
// FUNCIÓN PARA CALCULAR UNA INTEGRAL POR EL MÉTODO DEL RECTÁNGULO
// ###############################################################
real mn_rectangulo(
real (*f)(const real x), // Función que se integral
const real a, // Extremo izquierdo del intervalo
const real b, // Extremo derecho del intervalo
const int N) // Número de intervalos para calcular la integral
{
//HACER ALUMNO
real integral = 0; // Variable donde se acumula el valor de la integral
real funcion = (*f)(a); // Evaluación de la función en a
real h = (b-a); // Factor multiplicativo
for (int i = 0; i < N; i++) { // Bucle para acumular el valor de la integral
real ext_izq = a + i * (b-a)/N; // Extremo izquierdo subintervalo
real ext_dch = a + (i+1) * (b-a)/N; // Extremo derecho subintervalo
real ev_funcion = (*f)(ext_dch); // Evaluación de la función en el extremo derecho
integral += h * ((*f)((ext_izq + ext_dch)/2.)); // Acumulación valor integral
funcion = ev_funcion; // Actualización valor integral extremo izquierdo
}
return integral; // Se devuelve el valor de la integral
}
/************************************************************
FUNCION PARA CALCULAR UNA INTEGRAL POR EL METODO DE CUADRATURA DE GAUSS UTILIZANDO
UN VECTOR DE PUNTOS Y PESOS PARA EL INTERVALO [-1,1]. LOS VALORES DE PUNTOS Y PESOS
VIENEN EN LAS MATRICES w Y x, POR EJEMPLO PARA LA FÓRMULA DE 2 PUNTOS LOS PESOS
SERÍAN w[1][0] Y w[1][1] y LOS PUNTOS SERÍAN x[1][0] Y x[1][1]
NOTA : HAY QUE TENER EN CUENTA EL CAMBIO DE VARIABLE VISTO EN CLASE PARA PASAR UNA
INTEGRAL EN UN INTERVALO [a,b] AL INTERVALO [-1,1]
EN ESTE CASO LA INTEGRAL ES EN 2 DIMENSIONES Y POR TANTO LA FUNCIÓN A INTEGRAR DEPENDE
DE 2 VARIABLES Y HAY QUE AJUSTAR LA INTEGRAL TAL Y COMO SE EXPLICÓ EN CLASE. EL Nº DE
PUNTOS DE CUADRATURA A UTILIZAR ES EL MISMO EN LA VARIABLE x E y Y VIENE DADO EL PARÁMETRO
DE ENTRADA N
************************************************************/
real mn_cuadratura(
real (*f)(real x,real y) /* funcion que se integra */,
real a /* extremo izquierdo del intervalo para x*/,
real b /* extremo derecho del invervalo para x*/,
real c /* extremo inferior del intervalo para y*/,
real d /* extremo superior del invervalo para y*/,
Array2D< real > &x, /* matriz con los puntos de la tabla de cuadratura para el intervalo [-1,1]*/
Array2D< real > &w, /* matriz con los pesos de la tabla de cuadratura para el intervalo [-1,1] */
int N) /* Nº de puntos a usar en la fórmula de cuadratura */
{
// HACER ALUMNO
real suma = 0.; // variable donde guardare el resultado
real paso = 0.; // variables auxiliares que creo para no tener una linea demasiado larga
real paso2 = 0.; // a la hora de hacer los calculos
real imagen = 0.;
for(int k = 0; k<N; k++){ // bucle de la primera integral
for(int j = 0; j<N; j++){ //bucle de la segunda integral
paso = ((b-a)*x[N-1][k]+b+a)/2.; // valor que depende del primer intervalo [a,b]
paso2 = ((d-c)*x[N-1][j]+d+c)/2.; // valor que depende del segundo intervalo [c,d]
imagen = (*f)(paso,paso2); // calculo la funcion en esos puntos
suma += w[N-1][k]*w[N-1][j]*((b-a)/2.)*((d-c)/2.)*imagen; // calculo el sumatorio
}
}
return suma;
}
/************************************************************
FUNCION PARA CALCULAR UNA INTEGRAL POR EL METODO DE CUADRATURA DE GAUSS UTILIZANDO
UN VECTOR DE PUNTOS Y PESOS PARA EL INTERVALO [-1,1]
NOTA : HAY QUE TENER EN CUENTA EL CAMBIO DE VARIABLE VISTO EN CLASE PARA PASAR UNA
INTEGRAL EN UN INTERVALO [a,b] AL INTERVALO [-1,1]
LA FUNCIÓN DEVUELVE EL VALOR DE LA INTEGRAL.
************************************************************/
real mn_cuadratura(
real (*f)(real x) /* funcion que se integra */,
real a /* extremo izquierdo del intervalo */,
real b /* extremo derecho del invervalo*/,
Array1D< real > &x, /* vector con los puntos de la tabla de cuadratura para el intervalo [-1,1]*/
Array1D< real > &w) /* vector con los pesos de la tabla de cuadratura para el intervalo [-1,1] */
{
//A IMPLEMENTAR POR EL ALUMNO.
real suma = 0.;
real imagen = 0.;
real paso = 0.;
int N = x.dim();
for (int k = 0; k<N; k++) {
paso = ((b-a)*x[k]+b+a)/2.;
imagen = (*f)(paso);
suma = suma + w[k]*((b-a)/2.)*imagen;
}
return suma;
}
//Funcion que realiza la aproximacion de una derivada con 3 ptos
real mn_derivada3(real (*f)(real x),real xr, real xi){
real h = xr-xi;
real derivada;
derivada =(f(xi+h)-f(xi-h))/2;
return derivada;
}
// Funcion que realiza la derivada con 2 ptos.
real mn_deriva2(real (*f)(real x),real xr, real xi){
if (xr == xi) {return -1.;}
real derivada = (f(xr)-f(xi))/(xr-xi);
return derivada;
}
//============================================================================================
// FUNCION QUE CALCULA UNA INTEGRAL POR ITERACIONES DEL METODO DE SIMPSON:
// - SE PARTE DE N=1 INTERVALOS Y SE CALCULA LA INTEGRAL CON LA FUNCION ANTERIOR (mn_simpson)
// - SE ACTUALIZA N MULTIPLICANDOLO POR 2 Y SE CALCULA DE NUEVO LA INTEGRAL
// - SE REALIZA ESTE PROCESO HASTA QUE LA DISTANCIA RELATIVA ENTRE LOS 2 VALORES DE LA INTEGRAL
// SEA MENOR QUE LA TOLERANCIA O SE ALCANCE EL NÚMERO MÁXIMO DE INTERVALOS.
// - EN Niter SE ALMACENA EL NUMERO DE ITERACIONES QUE SE REALIZAN EN ESTE PROCESO.
//============================================================================================
real mn_simpson(
real (*f)(real x) /* funcion que se integra */,
real a /* extremo izquierdo del intervalo */,
real b /* extremo derecho del invervalo*/,
real TOL /* tolerancia para controlar la convergencia */,
real Imax /* número máximo de intervalos */,
int &Niter /* variable de salida con el numero de iteraciones realizadas */
)
{
/* A IMPLEMENTAR POR EL ALUMNO */
//NOTA: CALCULAR EL ERROR CON mn_distancia()
real integral = 0.;
real integral1 = 0.;
Niter=0;
real x= (b-a)/Imax;
for (int j=1;j<Imax;j++){
integral= mn_simpson(f,a,b,j);
integral1= mn_simpson(f,a,b,(2*j));
Niter++;
if (mn_distancia(integral,integral1)<TOL || j==Imax-1){
break;
}
}
}