Czy można zoptymalizować ten kod integracji, aby działał szybciej?

double trap(double func(double), double b, double a, double N) {
  double j;
  double s;
  double h = (b-a)/(N-1.0); //Width of trapezia

  double func1 = func(a);
  double func2;

  for (s=0,j=a;j<b;j+=h){
    func2 = func(j+h);
    s = s + 0.5*(func1+func2)*h;
    func1 = func2;
  }

  return s;
}

Powyżej jest mój kod C ++ dla integracji numerycznej 1D (przy użyciu reguły rozszerzonego trapezu) func()między granicami$[a,b]$ za pomocą $N-1$ trapezia.

W rzeczywistości wykonuję integrację 3D, w której ten kod jest wywoływany rekurencyjnie. Ja pracuję z$N = 50$ dając mi przyzwoite wyniki.

Inne niż redukcja $N$ponadto, czy ktoś jest w stanie zasugerować, jak zoptymalizować powyższy kod, aby działał szybciej? A może nawet sugeruje szybszą metodę integracji?

Matematycznie twoje wyrażenie jest równoważne z:

Abyś mógł to zaimplementować. Jak już powiedziano, czas prawdopodobnie jest zdominowany przez ocenę funkcji, więc aby uzyskać tę samą dokładność, można zastosować lepszą metodę integracji, która wymaga mniejszej oceny funkcji.

Kwadratura Gaussa jest we współczesnych czasach czymś więcej niż zabawką; przydatne tylko, jeśli potrzebujesz bardzo niewielu ocen. Jeśli chcesz czegoś łatwego do wdrożenia, możesz użyć reguły Simpsona, ale nie poszedłbym dalej niż zamówienie bez dobrego powodu.$1/N^3$

Jeśli krzywizna funkcji bardzo się zmienia, możesz użyć procedury kroku adaptacyjnego, która wybrałaby większy krok, gdy funkcja jest płaska, i mniejszy, bardziej dokładny, gdy krzywizna jest wyższa.

Możliwe, że ocena funkcji jest najbardziej czasochłonną częścią tego obliczenia. W takim przypadku powinieneś skupić się na poprawie prędkości func (), zamiast próbować przyspieszyć samą procedurę integracji.

W zależności od właściwości func () jest również prawdopodobne, że można uzyskać bardziej precyzyjną ocenę całki przy mniejszej liczbie ocen funkcji za pomocą bardziej wyrafinowanej formuły całkowania.

Możliwy? Tak. Przydatny? Nie. Optymalizacje, które wymienię tutaj, raczej nie spowodują więcej niż ułamek procenta różnicy w czasie wykonywania. Dobry kompilator może już to dla Ciebie zrobić.

W każdym razie, patrząc na twoją wewnętrzną pętlę:

    for (s=0,j=a;j<b;j+=h){
        func2 = func(j+h);
        s = s + 0.5*(func1+func2)*h;
        func1 = func2;
    }

Przy każdej iteracji pętli wykonujesz trzy operacje matematyczne, które można wyprowadzić na zewnątrz: dodawanie j + h, mnożenie przez 0.5i mnożenie przez h. Pierwszą, którą możesz naprawić, uruchamiając zmienną iteratora a + h, a drugą, dzieląc mnożenie:

    for (s=0, j=a+h; j<=b; j+=h){
        func2 = func(j);
        s += func1+func2;
        func1 = func2;
    }
    s *= 0.5 * h;

Chociaż chciałbym zauważyć, że robiąc to, z powodu błędu zaokrąglenia zmiennoprzecinkowego można pominąć ostatnią iterację pętli. (Było to również problem w swojej pierwotnej implementacji). Aby obejść to, używać unsigned intlub size_tlicznik:

    size_t n;
    for (s=0, n=0, j=a+h; n<N; n++, j+=h){
        func2 = func(j);
        s += func1+func2;
        func1 = func2;
    }
    s *= 0.5 * h;

Jak mówi odpowiedź Briana, lepiej poświęcić swój czas na optymalizację oceny funkcji func. Jeśli dokładność tej metody jest wystarczająca, wątpię, aby znaleźć coś szybciej dla tego samego N. (Chociaż możesz przeprowadzić kilka testów, aby sprawdzić, czy np. Runge-Kutta pozwala obniżyć się na Ntyle, że ogólna integracja zajmuje mniej czasu bez poświęcania dokładności.)

Jest kilka zmian, które poleciłbym ulepszyć obliczenia:

• Aby uzyskać wydajność i dokładność, użyj funkcji std::fma(), która wykonuje stopione dodawanie wielokrotne .
• Aby zwiększyć wydajność, należy odłożyć pomnożenie powierzchni każdego trapezu o 0,5 - możesz to zrobić raz na końcu.
• Unikaj wielokrotnego dodawania h, ponieważ może to powodować błędy zaokrąglania.

Ponadto wprowadziłbym kilka zmian dla przejrzystości:

• Nadaj tej funkcji bardziej opisową nazwę.
• Zamień kolejność ai bw podpisie funkcji.
• Zmień nazwę N→ n, h→ dx, j→ x2, s→ accumulator.
• Zmień nna int.
• Deklaruj zmienne w bardziej ścisłym zakresie.

#include <cmath>

double trapezoidal_integration(double func(double), double a, double b, int n) {
    double dx = (b - a) / (n - 1);   // Width of trapezoids

    double func_x1 = func(a);
    double accumulator = 0;

    for (int i = 1; i <= n; i++) {
        double x2 = a + i * dx;      // Avoid repeated floating-point addition
        double func_x2 = func(x2);
        accumulator = std::fma(func_x1 + func_x2, dx, accumulator); // Fused multiply-add
        func_x1 = func_x2;
    }

    return 0.5 * accumulator;
}

Jeśli twoja funkcja jest wielomianem, prawdopodobnie ważonym przez jakąś funkcję (np. Gaussa), możesz wykonać dokładną integrację w 3d bezpośrednio z formułą kubatury (np. Http://people.sc.fsu.edu/~jburkardt/c_src/ stroud / stroud.html ) lub z rzadką siatką (np. http://tasmanian.ornl.gov/ ). Metody te po prostu określają zestaw punktów i wag, przez które należy pomnożyć wartość funkcji, dzięki czemu są one bardzo szybkie. Jeśli twoja funkcja jest wystarczająco płynna, aby przybliżać ją wielomianami, wówczas metody te mogą dać bardzo dobrą odpowiedź. Formuły specjalizują się w typie funkcji, którą integrujesz, więc znalezienie odpowiedniego może wymagać trochę czasu.

Kiedy próbujesz obliczyć całkę numerycznie, próbujesz uzyskać pożądaną precyzję przy jak najmniejszym wysiłku, lub alternatywnie, próbujesz uzyskać najwyższą możliwą precyzję przy stałym wysiłku. Wygląda na to, że pytasz, jak sprawić, by kod dla jednego konkretnego algorytmu działał tak szybko, jak to możliwe.

To może dać ci mały zysk, ale będzie niewiele. Istnieją znacznie bardziej wydajne metody integracji numerycznej. Google za „Regułę Simpsona”, „Runge-Kutta” i „Fehlberg”. Wszystkie działają dość podobnie, oceniając niektóre wartości funkcji i sprytnie dodając wielokrotności tych wartości, wytwarzając znacznie mniejsze błędy przy tej samej liczbie ocen funkcji lub ten sam błąd przy znacznie mniejszej liczbie ocen.

Istnieje wiele sposobów integracji, z których reguła trapezowa jest najprostsza.

Jeśli wiesz cokolwiek na temat faktycznej funkcji, którą integrujesz, możesz to zrobić lepiej, jeśli ją wykorzystasz. Chodzi o to, aby zminimalizować liczbę punktów siatki w granicach dopuszczalnego poziomu błędu.

Na przykład trapezoidalnie dopasowuje się liniowo do kolejnych punktów. Możesz wykonać kwadratowe dopasowanie, które przy gładkiej krzywej pasowałoby lepiej, co pozwoliłoby na użycie grubszej siatki.

Symulacje orbit są czasem wykonywane przy użyciu stożków, ponieważ orbity są bardzo podobne do przekrojów stożkowych.

W mojej pracy integrujemy kształty zbliżone do krzywych w kształcie dzwonu, więc skuteczne jest ich modelowanie w ten sposób ( adaptacyjna kwadratura gaussowska jest w tej pracy uważana za „złoty standard”).

Tak więc, jak wskazano w innych odpowiedziach, zależy to w dużej mierze od tego, jak droga jest twoja funkcja. Optymalizacja kodu trapz jest tego warta tylko wtedy, gdy jest to naprawdę wąskie gardło. Jeśli nie jest to całkowicie oczywiste, powinieneś to sprawdzić, profilując swój kod (narzędzia takie jak Intels V-tune, Valgrind lub Visual Studio mogą to zrobić).

Sugerowałbym jednak zupełnie inne podejście: integracja z Monte Carlo . Tutaj po prostu przybliżasz całkę, próbkując swoją funkcję w losowych punktach, dodając wyniki. Zobacz ten pdf oprócz strony wiki, aby uzyskać szczegółowe informacje.

Działa to wyjątkowo dobrze w przypadku danych o dużych wymiarach, zwykle znacznie lepszych niż metody kwadraturowe stosowane w integracji 1-d.

Prosty przypadek jest bardzo łatwy do wdrożenia (patrz pdf), należy tylko uważać, aby standardowa funkcja losowa w c ++ 98 była dość zła zarówno pod względem wydajności, jak i jakości. W c ++ 11 możesz użyć Mersenne Twister w.

Jeśli twoja funkcja ma wiele odmian w niektórych obszarach, a mniej w innych, rozważ zastosowanie próbkowania warstwowego. Radziłbym raczej korzystać z biblioteki naukowej GNU , niż pisać własną.

W rzeczywistości wykonuję integrację 3D, w której ten kod jest wywoływany rekurencyjnie.

kluczem jest „rekurencyjnie”. Albo przechodzisz przez duży zestaw danych i rozważasz wiele danych więcej niż jeden raz, albo faktycznie generujesz swój zestaw danych na podstawie funkcji (częściowych?).

Rekursywnie oceniane integracje będą absurdalnie drogie i absurdalnie nieprecyzyjne, gdy moce wzrośnie w rekurencji.

Utwórz model interpolacji zestawu danych i dokonaj częściowej integracji symbolicznej. Ponieważ wiele danych ulega następnie zawaleniu na współczynniki funkcji bazowych, złożoność głębszej rekurencji rośnie bardziej wielomianowo (a zwykle raczej słabo), a nie wykładniczo. Otrzymujesz „dokładne” wyniki (nadal musisz opracować dobre schematy oceny, aby uzyskać rozsądną wydajność numeryczną, ale nadal powinno być wykonalne uzyskanie lepszej niż integracji trapezoidalnej).

Jeśli spojrzysz na oszacowania błędów dla reguł trapezoidalnych, zobaczysz, że są one powiązane z pewną pochodną zaangażowanych funkcji, a jeśli integracja / definicja zostanie wykonana rekurencyjnie, funkcje nie będą miały dobrze zachowanych pochodnych .

Jeśli twoim jedynym narzędziem jest młotek, każdy problem wygląda jak gwóźdź. Chociaż ledwo dotykasz problemu w swoim opisie, podejrzewam, że rekurencyjne stosowanie reguły trapezoidalnej jest niewłaściwe: dostajesz eksplozję zarówno niedokładności, jak i wymagań obliczeniowych.

oryginalny kod ocenia funkcję w każdym z N punktów, następnie dodaje wartości i mnoży sumę przez rozmiar kroku. jedyną sztuczką jest to, że wartości na początku i na końcu są dodawane z ciężarem , podczas gdy wszystkie punkty w środku są dodawane z pełnym ciężarem. w rzeczywistości są one również dodawane o wadze ale dwukrotnie. zamiast dodawać je dwa razy, dodaj je tylko raz z pełną wagą. oblicz mnożenie przez wielkość kroku poza pętlą. to wszystko, co można zrobić, aby to przyspieszyć, naprawdę.$1/2$$1/2$

    double trap(double func(double), double b, double a, double N){
double j, s;
double h = (b-a)/(N-1.0); //Width of trapezia

double s = 0;
j = a;
for(i=1; i<N-1; i++){
  j += h;
  s += func(j);
}
s += (func(a)+func(b))/2;

return s*h;
}