Cómo usar Matlab para evaluar [matemáticas] \ int ^ 1_0 \ int ^ 1_0 \ int ^ 1_0 \ frac {1} {(xyz) ^ {xyz}} ~ \ mathrm {d} x ~ \ mathrm {d} y ~ \ mathrm {d} z [/ math]

¡Hay muchas opciones! Aquí hay algunos:

Integración numérica ingenua

% De áreas de rectángulos:
[correo electrónico protegido] (x, y, z) (x. * y. * z). ^ (- x. * y. * z); % de la función que desea integrar
dx = 1e-2; % paso de discretización
[X, Y, Z] = ndgrid (dx / 2: dx: 1); % Omito los puntos finales izquierdos para evitar 0 ^ 0 issuses
F = f (X, Y, Z);
my_integral = sum (F (:)) * dx ^ 3; % la respuesta es aproximadamente 1.215

El problema con este enfoque numérico es que se encuentra con problemas de memoria muy rápidamente. Por ejemplo, el requisito de almacenamiento de lo que escribí es proporcional a dx ^ (- 3) = 10 ^ 6. Pero si hubiera hecho dx = 1e-3, habría necesitado 1000 veces el almacenamiento que excedería las limitaciones de memoria. Además, si queremos una integral cuádruple en lugar de una integral triple, el almacenamiento se escala como dx ^ -4 = 10 ^ 8 para dx = 1e-2. Eso está en el límite superior de lo que se puede hacer en mi máquina.


Monte Carlo

[correo electrónico protegido] (x, y, z) (x. * y. * z). ^ (- x. * y. * z); % de la función que desea integrar
n_trials = 1e7; % número de puntos iid para muestrear
r = rand (n_trials, 3); % de muestras de iid uniformes en un cubo
my_integral2 = mean (f (r (:, 1), r (:, 2), r (:, 3))); % la respuesta es aproximadamente 1.215

Este enfoque es bueno porque se escala muy bien con la dimensión. Si su integral hubiera sido de 4 dimensiones, podríamos usar el mismo número de puntos y aún así obtener una buena aproximación.


Simbólico
Este enfoque requiere la caja de herramientas simbólicas y todavía falla en este problema, pero lo escribiré de todos modos ya que funcionaría para problemas que el integrador simbólico puede entender.

syms xyz
my_integral3 = int (int (int ((x * y * z) ^ (- x * y * z), x, 0,1), y, 0,1), z, 0,1); % falla


Triplequad incorporado

[correo electrónico protegido] (x, y, z) (x. * y. * z). ^ (- x. * y. * z); % de la función que desea integrar
tol = 1e-6; %tolerancia
my_integral4 = triplequad (f, 0,1,0,1,0,1, tol)% devuelve 1.214832


Probablemente hay otras formas, pero esto debería ayudarlo a comenzar.

Sugiero instruir a Matlab para que calcule una suma aproximada de la serie

[matemáticas] \ begin {align} \ frac {1} {2} \ sum_ {k \ geq 0} \ frac {k ^ 2 + 3k + 2} {(k + 1) ^ {k + 3}}, \ end {align} [/ math]

porque la suma de esta serie es el valor de tu integral. Derivaré esto, a través de la expansión de serie exponencial estándar dentro de la integral como lo sugiere Yan Robert (por lo que mi contribución al hilo es solo para mostrar que la serie en realidad se reduce a la forma compacta y bonita escrita anteriormente).

Como escribió Yan, la función [matemática] 1 / (xyz) ^ {xyz} [/ matemática] se da de forma puntual en el cubo por la suma de la serie

[matemáticas] \ begin {align} & \ frac {1} {(xyz) ^ {xyz}} \\ & \ quad \ quad = e ^ {- xyz \ log (xyz)} \\ & \ quad \ quad = \ sum_ {k \ geq 0} \ frac {(- xyz \ log (xyz)) ^ k} {k!}. \ end {align} [/ math]

Esta es una suma de términos positivos (ya que [math] xyz <1 [/ math] en el cubo), entonces, por convergencia monótona, la integral de la suma es igual a la suma de las integrales:

[matemáticas] \ begin {align} & \ int _ {(0,1) ^ 3} \ frac {1} {(xyz) ^ {xyz}} dxdydz \\ & = \ sum_ {k \ geq 0} \ frac { (-1) ^ k} {k!} \ Underbrace {\ int _ {(0,1) ^ 3} (xyz) ^ k (\ log xyz) ^ kdxdydz} _ {I (k)}. \ end {align} [/ math]

Ahora, para calcular la integral [matemáticas] I (k) = \ int _ {(0,1) ^ 3} (xyz) ^ k (\ log xyz) ^ kdxdydz [/ math] argumentamos lo siguiente. Si [math] \ lambda [/ math] es un parámetro complejo, defina

[matemáticas] \ begin {align} f (\ lambda) = \ int _ {(0,1) ^ 3} (xyz) ^ \ lambda dxdydz = (\ lambda + 1) ^ {- 3}. \ end {align} [/ math]

Evidentemente, si [math] f ^ {(j)} [/ math] denota la [math] j [/ math] -th derivada de [math] f [/ math] entonces

[matemáticas] \ begin {align} \ int _ {(0,1) ^ 3} (xyz) ^ \ lambda (\ log xyz) ^ jdxdydz = f ^ {(j)} (\ lambda). \ end {align} [/matemáticas]

La diferenciación dentro de la integral es bastante fácil de justificar en este caso ya que el logaritmo tiene una singularidad tan débil en cero. Alternativamente, como solo estamos interesados ​​en valores enteros no negativos de [math] \ lambda [/ math], los valores de las integrales correspondientes se pueden calcular directamente usando la integración por partes.

De todos modos, sabemos que [math] f (\ lambda) = (\ lambda + 1) ^ {- 3} [/ math] para que podamos diferenciar esta función racional directamente para obtener

[matemáticas] \ begin {align} f ^ {(j)} (\ lambda) = \ frac {(- 1) ^ j (j + 2)!} {2! (\ lambda + 1) ^ {j + 3 }}. \ end {align} [/ math]

Por lo tanto, [matemáticas] I (k) = f ^ {(k)} (k) [/ matemáticas] y

[matemáticas] \ begin {align} & \ int _ {(0,1) ^ 3} \ frac {dxdydz} {(xyz) ^ {xyz}} \\ & = \ sum_ {k \ geq 0} \ frac {( -1) ^ k} {k!} I (k) \\ & = \ sum_ {k \ geq 0} \ frac {(- 1) ^ k} {k!} F ^ {(k)} (k) \\ & = \ sum_ {k \ geq 0} \ frac {(- 1) ^ k} {k!} \ frac {(- 1) ^ k (k + 2)!} {2! (k + 1) ^ {k + 3}} \\ & = \ frac {1} {2} \ sum_ {k \ geq 0} \ frac {k ^ 2 + 3k + 2} {(k + 1) ^ {k + 3} } \ end {align} [/ math]

como se desee.

Incidentalmente, para la integral dimensional [matemática] d [/ matemática] análoga, obtendríamos fácilmente

[matemáticas] \ begin {align} \ int _ {(0,1) ^ d} & (x_1 \ cdots x_d) ^ \ lambda (\ log x_1 \ cdots x_d) ^ jdx_1 \ cdots dx_d \\ & = \ left (\ frac {d} {d \ lambda} \ right) ^ j \ frac {1} {(\ lambda + 1) ^ d} \\ & = \ frac {(- 1) ^ j (j + d-1)! } {(d-1)! (\ lambda + 1) ^ {j + d}} \ end {align} [/ math]

y por lo tanto

[matemáticas] \ begin {align} & I (k, d) \\ & = \ int _ {(0,1) ^ d} (x_1 \ cdots x_d) ^ k (\ log x_1 \ cdots x_d) ^ kdx_1 \ cdots dx_d \\ & = \ frac {(- 1) ^ k (k + d-1)!} {(d-1)! (k + 1) ^ {k + d}} \ end {align} [/ math]

así que eso

[matemáticas] \ begin {align} & \ int _ {(0,1) ^ d} \ frac {dx_1 \ cdots dx_d} {(x_1 \ cdots x_d) ^ {x_1 \ cdots x_d}} \\ & = \ sum_ { k \ geq 0} \ frac {(- 1) ^ k} {k!} I (k, d) \\ & = \ sum_ {k \ geq 0} \ frac {(- 1) ^ k} {k! } \ frac {(- 1) ^ k (k + d-1)!} {(d-1)! (k + 1) ^ {k + d}} \\ & = \ frac {1} {(d -1)!} \ Sum_ {k \ geq 0} \ frac {(k + d-1)!} {K! (K + 1) ^ {k + d}}. \ end {align} [/ math]

Esta es una integral interesante … ¡gracias por publicar!

Enfoque numérico: ¿cómo muestrear?
Aquí se han sugerido dos métodos para el enfoque numérico, muestreo regular y muestreo aleatorio.
Sin embargo, en muchos casos, el muestreo semialeatorio tiene una convergencia superior.

El problema del muestreo de cuadrícula regular es que muestreas en frecuencias espaciales fijas, con el riesgo de submuestreo / sobremuestreo de frecuencias importantes.

El problema del muestreo aleatorio es que podría sub / muestrear fácilmente regiones importantes por casualidad.

Un mejor enfoque es el muestreo semialeatorio , en el que los puntos muestreados se distribuyen globalmente de manera más uniforme.
En este ejemplo, elegí muestrear cerca de cada punto en la cuadrícula regular, sin embargo, la distancia relativa a este punto de la cuadrícula fue variada. En cada dimensión, seleccioné al azar la desviación de una distribución uniforme en el intervalo [-1 / 2 * d , 1/2 * d ] en el que d es la distancia entre los puntos de la cuadrícula regular.

En la siguiente figura se ve cómo el muestreo semialeatorio (línea azul) converge más rápido que el muestreo aleatorio (línea roja) o el muestreo regular (línea verde).

El muestreo semi aleatorio (línea azul) converge más rápido que el muestreo aleatorio (línea roja) o el muestreo regular (línea verde). El eje horizontal (N) representa el número de puntos de cuadrícula en cada dimensión, lo que lleva a N ^ 3 puntos muestreados. (También N ^ 3 muestras para muestreo aleatorio).

La misma cifra, pero ahora para diez simulaciones.

%% CODIGO MATLAB:
Nseries = 50: 50: 400;
resultados = ceros (3, longitud (Nseries));
para m = 1: longitud (Nseries)

N = Nseries (m);

% Muestreo aleatorio
xyz = rand (1, N ^ 3). * rand (1, N ^ 3). * rand (1, N ^ 3);
resultados (1, m) = media (1./xyz.^xyz);

% De muestreo regular
gridPoints_1D = linspace (1 / (2 * N), 1-1 / (2 * N), N);
[x, y, z] = meshgrid (gridPoints_1D);

% en nodos de cuadrícula exactos:
xyz = x (:). * y (:). * z (:);
resultados (2, m) = media (1./xyz.^xyz);

% cerca de nodos de cuadrícula, desviación aleatoria:
x = x + (rand (tamaño (x)) – 0.5) / N;
y = y + (rand (tamaño (x)) – 0.5) / N;
z = z + (rand (tamaño (x)) – 0.5) / N;
xyz = x (:). * y (:). * z (:);
resultados (3, m) = media (1./xyz.^xyz);
fin

figura;
semilogx (Nseries, resultados (1, :), ‘r’); Espere;
semilogx (Nseries, resultados (2, :), ‘g’);
semilogx (Nseries, resultados (3, :), ‘b’);

Para otros esquemas de muestreo semi aleatorio, ver:
Método cuasi-Monte Carlo

MATLAB tiene la opción de calcular integrales triples. Evaluar numéricamente la integral triple

Crear una función

function f = tripinteg(x,y,z) f = 1./((x*y*z).^(x*y*z)); end

function f = tripinteg(x,y,z) f = 1./((x*y*z).^(x*y*z)); end

function f = tripinteg(x,y,z) f = 1./((x*y*z).^(x*y*z)); end

Y luego corre

Q = triplequad(@tripinteg,0,1,0,1,0,1);

Alternativamente, también podría buscar la respuesta de Justin Rising

Monte Carlo lo. Tome una gran cantidad de muestras de la distribución uniforme en [matemáticas] [0, 1] [/ matemáticas], conéctelas a la función y simplemente promedie los resultados. Esto le dará una muy buena aproximación a la integral.

Use la función integral3 en matlab para evaluar integrales triples

fun = @ (x, y, z) 1./((x.*y.*z).^(x.*y.*z))% Define la función anónima para la ecuación dada%
% Calcular la integral triple con los límites provistos
q = integral3 (diversión, 0,1,0,1,0,1)
%responder
q = 1.21483799428875

Usaría series de Taylor para escribir

[matemáticas] \ int_0 ^ 1 \ int_0 ^ 1 \ int_0 ^ 1 e ^ {- xyz \ log (xyz)} dx dy dz = \ int_0 ^ 1 \ int_0 ^ 1 \ int_0 ^ 1 \ sum_ {n = 0} ^ \ infty (-xyz \ log (xyz)) ^ n / n !. [/matemáticas]

Como xyz <= 1 significa que [math] xyz \ log (xyz) [\ math] está acotado y se aproxima a cero como lo hace xyz, la serie anterior convergerá. Puede truncar lo anterior con la precisión deseada y luego integrar manualmente la nueva expresión.