[C++] passare da integrazione unidimensionale a multidimensionale

[alastor1984]

Salve.
Per la mia tesi sto facendo dei calcoli in c++.
Ho scritto una funzione che, prendendo in pasto un'altra funzione e dei parametri, mi calcola l'integrale improprio tra -\inf e + \inf.
In linea di principio, dovrebbe essere facile adattare questo algoritmo al caso multidimensionale, dato che con fubini tonelli basta fissare N-1 variabili, integrarne 1 e ricorsivamente chiamare la funzione finchè non le saturi tutte..
Però non riesco a scrivere tutto questo in codice, qualcuno ha delle idee?
Il sistema operativo è Mac Osx, l'Ide è Eclipse e il compilatore G++.
ecco il codice che ho scritto.
Grazie mille,
Enzo

#ifndef INTEGRAL_
#define INTEGRAL_
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <math.h>
#include <stdio.h> 
#include "tnt.h"

using namespace std;
using namespace TNT;

double change(double (*fun)(double), double x)
{	double y=sinh(M_PI_2*sinh(x));
	double dydx=M_PI_2*cosh(M_PI_2*sinh(x))*cosh(x);
	return fun(y)*dydx;
	}; 

double integral1D(double (*fun)(double))
{//mode =0: trapezium rule
 //mode =1: simpson rule
	int mode=1;
	int NMAX=20;
	double a=-4.0, b=4.0;
	double h=(b-a);
	double toll=pow(10, -6.0);
	double funca=change(fun,a),funcb=change(fun,b);
	double onethird=1.0/3.0;
	TNT::Array1D<double> sum(NMAX,0.0);
	TNT::Array1D<double> estimate(NMAX,0.0);
	sum[0]=(funca+funcb)*0.5;
	estimate[0]=sum[0];
	// cout << "sum="<<  sum[0] <<endl;
	for(int n=1; n<NMAX; n++)
	{
		double x=a+h*0.5;
		//cout << "h="<< h <<endl;
		// cout << "x=" << x<< endl;
		for (;x<b; x+=h){
			
			sum[n]+=change(fun,x);
			//cout <<"x= "<< x <<  " func(x)= " <<func(x) << endl;
		}
		sum[n]=sum[n]+sum[n-1];
		h=h*0.5;
		if(mode==0) 
					{
					estimate[n]=(4.0*sum[n]-sum[n-1]*2)*h*onethird;
					//cout << "estimate at step "<< n <<"="<<estimate[n]<< endl;
					if (n>5 && abs(estimate[n]-estimate[n-1])< toll)	
					return estimate[n];
					}
		if(mode==1) 
					{
					estimate[n]=sum[n]*h;
					//cout << "estimate at step "<< n <<"="<<estimate[n]<< endl;
				
					if (n>5 && abs(estimate[n]-estimate[n-1])< toll)
				
					return estimate[n];
					}
		//std::cout << "prova= " << sum << std::endl; 
}cout <<  " the desired precision cannot be reached with requested number of steps:";
};



#endif /*INTEGRAL_*/

[alastor1984]

si certo! dovrebbe venire 0 per ragioni di simmetria..
Per il tipo di applicazioni che uso (trasformate di p.d.f. ) l'integrale converge..
per funzioni che hanno "esplosioni in mezzo, devi usare diversi accorgimenti, che nel mio caso non sono adottati.

[alastor1984]

bhe veramente questa è integrazione di Riemann.. L'integrale di Lebesgue maneggia una infinità numerabile di discontinuità, e il problema della funzione che esplode c'è sempre, il risultato a cui ti riferisci è:
se |f(x)|<inf, mis(A)=0 => \int_(A)(f(x))dx=0;
se f(x) fa il botto non è detto che sia neanche finita..
Credo che gli integrali di Lebesgue non si possano fare nemmeno al calcolatore, sono solo degli strumenti teorici che facilitano il passaggio sotto gli operatori differenziali!
ti rimando al Kolmogorov Fomin per maggiori dettagli!
fine OT
Nessuno sa come si può automatizzare? sennò mi tocca fare montecarlo, che è lento come il cucco, e dato che uso questi integrali per valutare Runge Kutta Fehlberg può vuoler dire metterci le ere geologiche per ogni problema!