Apakah mungkin untuk mengoptimalkan kode integrasi ini agar berjalan lebih cepat?

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;
}

Di atas adalah kode C ++ saya untuk integrasi numerik 1D (menggunakan aturan trapezium diperpanjang) dari func()antara batas $[a,b]$ menggunakan $N-1$ trapezia.

Saya sebenarnya melakukan integrasi 3D, di mana kode ini disebut secara rekursif. Saya bekerja dengan $N = 50$ memberi saya hasil yang layak.

Selain mengurangi lebih lanjut, adakah yang bisa menyarankan cara mengoptimalkan kode di atas sehingga berjalan lebih cepat? Atau, bahkan, dapat menyarankan metode integrasi yang lebih cepat?$N$

Secara matematis, ekspresi Anda setara dengan:

Jadi Anda bisa menerapkannya. Seperti yang dikatakan, waktu mungkin didominasi oleh evaluasi fungsi, jadi untuk mendapatkan akurasi yang sama, Anda dapat menggunakan metode integrasi yang lebih baik yang membutuhkan evaluasi fungsi yang lebih sedikit.

Quadrature Gaussian, di zaman modern, lebih dari sekadar mainan; hanya berguna jika Anda memerlukan sangat sedikit evaluasi. Jika Anda menginginkan sesuatu yang mudah diimplementasikan, Anda dapat menggunakan aturan Simpson, tetapi saya tidak akan melangkah lebih jauh dari memesan tanpa alasan yang bagus.$1/N^3$

Jika kelengkungan fungsi banyak berubah, Anda bisa menggunakan langkah rutin adaptatif, yang akan memilih langkah yang lebih besar ketika fungsi rata, dan yang lebih kecil lebih akurat ketika kelengkungan lebih tinggi.

Kemungkinannya adalah bahwa evaluasi fungsi adalah bagian yang paling memakan waktu dari perhitungan ini. Jika itu masalahnya, maka Anda harus fokus pada peningkatan kecepatan func () daripada mencoba mempercepat rutinitas integrasi itu sendiri.

Bergantung pada properti func (), kemungkinan Anda juga bisa mendapatkan evaluasi integral yang lebih tepat dengan evaluasi fungsi yang lebih sedikit dengan menggunakan rumus integrasi yang lebih canggih.

Bisa jadi? Iya. Berguna? Tidak. Optimalisasi yang akan saya daftarkan di sini tidak mungkin menghasilkan lebih dari sepersekian kecil perbedaan persen dalam runtime. Kompiler yang baik mungkin sudah melakukan ini untuk Anda.

Lagi pula, melihat loop batin Anda:

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

Di setiap pengulangan Anda melakukan tiga operasi matematika yang bisa dibawa keluar: menambah j + h, mengalikan dengan 0.5, dan mengalikan dengan h. Yang pertama Anda bisa perbaiki dengan memulai variabel iterator Anda di a + h, dan yang lainnya dengan memperhitungkan keluar perkalian:

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

Meskipun saya akan menunjukkan bahwa dengan melakukan ini, karena kesalahan roundoff titik mengambang adalah mungkin untuk melewatkan iterasi terakhir dari loop. (Ini juga merupakan masalah dalam implementasi awal Anda.) Untuk menyiasatinya, gunakan unsigned intatau size_tkonter:

    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;

Seperti jawaban Brian katakan, waktu Anda lebih baik dihabiskan untuk mengoptimalkan evaluasi fungsi func. Jika akurasi metode ini cukup, saya ragu Anda akan menemukan sesuatu yang lebih cepat untuk hal yang sama N. (Meskipun Anda dapat menjalankan beberapa tes untuk melihat apakah misalnya Runge-Kutta memungkinkan Anda Ncukup rendah sehingga keseluruhan integrasi membutuhkan waktu lebih sedikit tanpa mengorbankan akurasi.)

Ada beberapa perubahan yang saya rekomendasikan untuk meningkatkan perhitungan:

• Untuk kinerja dan keakuratan, gunakan std::fma(), yang melakukan penambahan-pengganda yang menyatu .
• Untuk kinerja, tunda mengalikan luas masing-masing trapesium dengan 0,5 - Anda dapat melakukannya sekali di akhir.
• Hindari penambahan berulang h, yang bisa mengakumulasi kesalahan pembulatan.

Selain itu, saya akan membuat beberapa perubahan untuk kejelasan:

• Beri fungsi nama yang lebih deskriptif.
• Tukar urutan adan bdalam tanda tangan fungsi.
• Ganti nama N→ n, h→ dx, j→ x2, s→ accumulator.
• Ubah nke int.
• Deklarasikan variabel dalam lingkup yang lebih ketat.

#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;
}

Jika fungsi Anda adalah polinomial, mungkin dibobot oleh beberapa fungsi (misalnya gaussian), Anda dapat melakukan integrasi tepat dalam 3d langsung dengan rumus cubature (mis. Http://people.sc.fsu.edu/~jburkardt/c_src/ stroud / stroud.html ) atau dengan grid tipis (mis. http://tasmanian.ornl.gov/ ). Metode-metode ini hanya menentukan satu set poin dan bobot untuk mengalikan nilai fungsi dengan, sehingga mereka sangat cepat. Jika fungsi Anda cukup halus untuk diperkirakan oleh polinomial, maka metode ini masih dapat memberikan jawaban yang sangat baik. Rumus khusus untuk jenis fungsi yang Anda integrasikan, sehingga mungkin perlu beberapa penggalian untuk menemukan yang tepat.

Saat Anda mencoba menghitung integral secara numerik, Anda mencoba untuk mendapatkan ketelitian yang Anda inginkan dengan upaya sekecil mungkin, atau sebagai alternatif, cobalah untuk mendapatkan ketelitian setinggi mungkin dengan upaya tetap. Anda sepertinya bertanya bagaimana membuat kode untuk satu algoritma tertentu berjalan secepat mungkin.

Itu mungkin memberi Anda sedikit keuntungan, tetapi itu akan sedikit. Ada banyak metode yang lebih efisien untuk integrasi numerik. Google untuk "aturan Simpson", "Runge-Kutta", dan "Fehlberg". Mereka semua bekerja sangat mirip dengan mengevaluasi beberapa nilai fungsi dan secara pintar menambahkan beberapa nilai tersebut, menghasilkan kesalahan yang jauh lebih kecil dengan jumlah evaluasi fungsi yang sama, atau kesalahan yang sama dengan jumlah evaluasi yang jauh lebih kecil.

Ada banyak cara untuk melakukan integrasi, di mana aturan trapesium adalah yang paling sederhana.

Jika Anda tahu apa-apa tentang fungsi aktual yang Anda integrasikan, Anda dapat melakukan lebih baik jika Anda mengeksploitasinya. Idenya adalah untuk meminimalkan jumlah titik kisi dalam tingkat kesalahan yang dapat diterima.

Misalnya, trapesium membuat kecocokan linier ke poin berurutan. Anda bisa membuat kuadrat pas, yang jika kurva mulus akan cocok lebih baik, yang memungkinkan Anda untuk menggunakan kisi yang lebih kasar.

Simulasi orbital kadang-kadang dilakukan menggunakan kerucut, karena orbit sangat mirip dengan bagian kerucut.

Dalam karya saya, kami mengintegrasikan bentuk yang mendekati kurva berbentuk lonceng, sehingga sangat efektif untuk memodelkannya seperti itu ( quadrature Gaussian adaptif dianggap sebagai "standar emas" dalam karya ini).

Jadi, seperti yang telah ditunjukkan dalam jawaban lain, ini sangat tergantung pada seberapa mahal fungsi Anda. Mengoptimalkan kode trapz Anda hanya layak jika itu benar-benar menjadi hambatan Anda. Jika tidak sepenuhnya jelas, Anda harus memeriksanya dengan membuat profil kode Anda (alat seperti Intels V-tune, Valgrind atau Visual Studio dapat melakukan ini).

Namun saya akan menyarankan pendekatan yang sama sekali berbeda: integrasi Monte Carlo . Di sini Anda cukup memperkirakan integral dengan mengambil sampel fungsi Anda pada titik-titik acak menambahkan hasilnya. Lihat pdf ini di samping halaman wiki untuk detailnya.

Ini berfungsi sangat baik untuk data dimensi tinggi, biasanya jauh lebih baik daripada metode quadrature yang digunakan dalam integrasi 1-d.

Kasus sederhana sangat mudah diimplementasikan (lihat pdf), hanya hati-hati bahwa fungsi acak standar di c ++ 98 cukup buruk baik dari segi kinerja dan kualitas. Di c ++ 11, Anda dapat menggunakan Mersenne Twister di.

Jika fungsi Anda memiliki banyak variasi di beberapa area dan kurang di area lainnya, pertimbangkan untuk menggunakan pengambilan sampel bertingkat. Saya akan merekomendasikan menggunakan perpustakaan ilmiah GNU , daripada menulis sendiri.

Saya sebenarnya melakukan integrasi 3D, di mana kode ini disebut secara rekursif.

"secara rekursif" adalah kuncinya. Anda baik melalui kumpulan data besar dan mempertimbangkan banyak data lebih dari sekali, atau Anda benar-benar menghasilkan set data Anda sendiri dari fungsi (secara terpisah?).

Integrasi yang dievaluasi secara rekursif akan menjadi sangat mahal, dan sangat tidak tepat karena kekuatan meningkat dalam rekursi.

Buat model untuk menginterpolasi kumpulan data Anda dan lakukan integrasi simbolik piecewise. Karena banyak data kemudian runtuh menjadi koefisien fungsi-fungsi dasar, kompleksitas untuk rekursi yang lebih dalam tumbuh secara polinomi (dan biasanya agak rendah) daripada secara eksponensial. Dan Anda mendapatkan hasil "tepat" (Anda masih perlu memikirkan skema evaluasi yang baik untuk mendapatkan kinerja numerik yang masuk akal, tetapi harus tetap layak untuk mendapatkan yang lebih baik daripada integrasi trapesium).

Jika Anda melihat pada perkiraan kesalahan untuk aturan trapesium, Anda akan menemukan bahwa mereka terkait dengan beberapa turunan dari fungsi yang terlibat, dan jika integrasi / definisi dilakukan secara rekursif, fungsi tidak akan cenderung memiliki turunan yang berperilaku baik. .

Jika satu-satunya alat Anda adalah palu, setiap masalah terlihat seperti paku. Meskipun Anda hampir tidak menyentuh masalah dalam uraian Anda, saya memiliki kecurigaan bahwa menerapkan aturan trapesium secara rekursif adalah kecocokan yang buruk: Anda mendapatkan ledakan persyaratan ketidaktepatan dan komputasi.

kode asli mengevaluasi fungsi di setiap titik N, kemudian menambahkan nilai-nilai ke atas, dan mengalikan jumlah dengan ukuran langkah. satu-satunya trik adalah bahwa nilai-nilai di awal dan akhir ditambahkan dengan bobot , sementara semua poin di dalam ditambahkan dengan bobot penuh. sebenarnya, mereka juga ditambahkan dengan berat tetapi dua kali. alih-alih menambahkannya dua kali, tambahkan hanya sekali dengan bobot penuh. faktor keluar perkalian dengan ukuran langkah di luar loop. itu saja yang bisa dilakukan untuk mempercepat ini, sungguh.1 /$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;
}