Эффективный способ вычисления интегралов? - PullRequest
Новая
       8

Эффективный способ вычисления интегралов?

0 голосов
/ 17 мая 2018

Я знаю, как использовать Монте-Карло для вычисления интегралов, но мне было интересно, можно ли использовать правило трапеции в сочетании с numpy для эффективности, чтобы получить тот же интеграл, я не уверен, какой из них самый быстрый или возможно ли последнее?

например для интеграции e**-x**2 > y я мог бы использовать метод Монте-Карло, как это: 

import numpy as np
import matplotlib.pyplot as plt

X = np.random.rand(500000,2)
X[:,0] = X[:,0]*4-2
J = np.where(X[:,1] < np.exp(-X[:,0]**2))[0]
Xp = X[:2000]
Ip = [i for i in range(len(Xp)) if i in J]
Inp = [i for i in range(len(Xp)) if i not in J]
plt.plot(Xp[Ip,0],Xp[Ip,1], 'bd', Xp[Inp,0],Xp[Inp,1], 'rd')
plt.show()

И это очень легко вычислить: 

print len(J) / 500000.0 * 4

Что дает: 

1.767784

В этом случае это было легко, но если интервалы не указаны, как в [a,b] , n, и я хочу создать функцию, то я думаю, что приведенный выше метод не очень эффективен, по крайней мере, мне так кажется.

Итак, мой вопрос: могу ли я интегрировать непрерывную функцию, например, cos(x)/x для определенного интервала, например, [a,b], в функцию с правилом трапеции?

И это лучше, чем метод, который я использовал здесь?

Любой совет приветствуется.

Ответы [ 2 ]

0 голосов
/ 13 октября 2018

Вы также можете использовать приближение Римана.Ниже приведен код на Java 

package math;

import java.util.Optional;
import java.util.function.*;
import java.util.stream.IntStream;
import static java.lang.Math.*;

public class IntegralJava8
{
    public interface Riemann extends BiFunction<Function<Double, Double>, Integer, 
    BinaryOperator<Double>> { }

    public static void main(String args[])
    {
        int N=100000;
        Riemann s = (f, n) -> (a, b) -> 
        IntStream.range(0, n)
        .mapToDouble(i -> f.apply(a + i * ((b - a) / n)) * ((b - a) / n)).sum();
        Optional<Double> gaussIntegral = 
            Optional.of(s.apply(x -> exp(-pow(x, 2)), N).apply(-1000.0, 1000.0));
        gaussIntegral.ifPresent(System.out::println);
    }
}

. В вышеприведенном классе он вычислит гауссовский интеграл от -infinity до бесконечности, равный квадратному корню из PI (1.772)

0 голосов
/ 17 мая 2018

Просто используйте scipy.integrate.quad: 

from scipy import integrate
from np import inf
from math import exp, sqrt, pi

res, errEstimate = integrate.quad(lambda x: exp(-x**2), -inf, +inf)

print(res)       #output: 1.7724538509055159
print(sqrt(pi))  #output: 1.7724538509055159

Последняя строка просто проверяет, что вычисленный интеграл действительно является квадратным корнем от Pi (это Гауссовский интеграл ).

...