Matlab Integración numérica multivariable

[Lohengrin]

Hola en un problema necesito calcular una integral definida multivariable con intervalos de integración infinitos, es decir

$\displaystyle \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} f(x_1,x_2, \ldots , x_n)\,dx_1 \cdots dx_n$$$

Donde n puede tomar valores del orden 80-100.
No he sido capaz de encontrar una solución analítica ni de descomponer la función $f$ como producto de funciones unidimensionales en las variables $$$x_i$$$.

Nunca he tenido que hacer algo así por lo que parto de cero. Sabéis de alguna toolbox de Matlab o algún otro paquete de software con el que pueda resolver numéricamente la integral. El resultado de la integral se utiliza en un programa que realiza otros muchos calculos y que está implementado en matlab. que Muchas gracias

[Lohengrin]

Re: Integración numérica multivariable

[FONT=verdana]Pongo un script de ejemplo integrando una normal multivariante. He usado esta función por sencillez y porque sabemos la solución que es uno.[/FONT]

[FONT=verdana]clc[/FONT]
[FONT=verdana]close all[/FONT]
[FONT=verdana]clear all[/FONT]
[FONT=verdana]pmax=10;[/FONT]

[FONT=verdana]t=zeros(pmax,1);[/FONT]
[FONT=verdana]result=zeros(pmax,1);[/FONT]
[FONT=verdana]sigma=2;[/FONT]

[FONT=verdana]for iter=1 : pmax[/FONT]
[FONT=verdana] str = [' Iteration ', num2str(iter), ' of ',num2str(pmax)];[/FONT]
[FONT=verdana] disp(str)[/FONT]
[FONT=verdana] p=iter;[/FONT]
[FONT=verdana] S=(sigma^2)*eye(p);[/FONT]
[FONT=verdana] invS=inv(S);[/FONT]
[FONT=verdana] detS=det(S);[/FONT]

[FONT=verdana] %x=[x1;...; xp][/FONT]
[FONT=verdana] x=sym('x',[p,1]);[/FONT]

[FONT=verdana] % x'= complex conjugate transpose[/FONT]
[FONT=verdana] % x.'= transpose without conjugation[/FONT]
[FONT=verdana] xT=x.';[/FONT]
[FONT=verdana] f=exp(-0.5*xT*invS*x);[/FONT]
[FONT=verdana] fk=f;[/FONT]
[FONT=verdana] tstart=tic;[/FONT]
[FONT=verdana] for k=1 : p[/FONT]
[FONT=verdana] fk=int(fk,x(k),-inf,inf);[/FONT]
[FONT=verdana] end[/FONT]
[FONT=verdana] t(iter)=toc(tstart);[/FONT]
[FONT=verdana] aux=1/sqrt((2*pi)^p);[/FONT]
[FONT=verdana] aux=aux/sqrt(detS);[/FONT]
[FONT=verdana] aux=aux*fk;[/FONT]
[FONT=verdana] result(iter)=aux;[/FONT]
[FONT=verdana]end[/FONT]

[FONT=verdana]figure[/FONT]
[FONT=verdana]plot(result)[/FONT]
[FONT=verdana]title('Result (Ideally always 1)')[/FONT]
[FONT=verdana]figure[/FONT]
[FONT=verdana]plot(t)[/FONT]
[FONT=verdana]title('time elapsed')


[/FONT][FONT=verdana]El problema es que el tiempo de computo crece demasiado con la dimensión p. Si alguien me puede dar algún consejo para acelerar el computo lo agradeceré.[/FONT]