সিম্পসনের সূত্র দ্বারা সমাকলন#
আমরা একটি নির্দিষ্ট সমাকলনের মান গণনা করতে যাচ্ছি
$$\int_a ^ b f (x) dx$$এখানে বর্ণিত সমাধানটি থমাস সিম্পসন তাঁর ১৭৪৩ সালের একটি অভিসন্দর্ভে প্রকাশ করেছিলেন।
সিম্পসনের সূত্র#
ধরি $n$ কোনো স্বাভাবিক সংখ্যা। আমরা সমাকলনের ব্যবধান $[a, b]$ কে $2n$ টি সমান অংশে ভাগ করি:
$$x_i = a + i h, ~~ i = 0 \ldots 2n,$$$$h = \frac {b-a} {2n}.$$এখন আমরা প্রতিটি ব্যবধান $[x_ {2i-2}, x_ {2i}]$, $i = 1 \ldots n$ এ আলাদাভাবে সমাকলন গণনা করি, এবং তারপর সকল মান যোগ করি।
ধরি আমরা পরবর্তী ব্যবধান $[x_ {2i-2}, x_ {2i}], i = 1 \ldots n$ বিবেচনা করছি। এতে ফাংশন $f(x)$ কে ৩টি বিন্দু $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$ দিয়ে যাওয়া একটি পরাবৃত্ত $P(x)$ দ্বারা প্রতিস্থাপন করি। এরকম একটি পরাবৃত্ত সবসময় বিদ্যমান এবং অনন্য; এটি বিশ্লেষণাত্মকভাবে পাওয়া যায়। উদাহরণস্বরূপ আমরা ল্যাগ্রাঞ্জ বহুপদী ইন্টারপোলেশন ব্যবহার করে এটি নির্মাণ করতে পারি। এখন শুধু এই বহুপদীর সমাকলন করা বাকি। যদি আপনি একটি সাধারণ ফাংশন $f$ এর জন্য এটি করেন, তাহলে একটি অসাধারণ সরল রাশি পাওয়া যায়:
$$\int_{x_ {2i-2}} ^ {x_ {2i}} f (x) ~dx \approx \int_{x_ {2i-2}} ^ {x_ {2i}} P (x) ~dx = \left(f(x_{2i-2}) + 4f(x_{2i-1})+(f(x_{2i})\right)\frac {h} {3} $$সকল ব্যবধানে এই মানগুলো যোগ করলে আমরা চূড়ান্ত সিম্পসনের সূত্র পাই:
$$\int_a ^ b f (x) dx \approx \left(f (x_0) + 4 f (x_1) + 2 f (x_2) + 4f(x_3) + 2 f(x_4) + \ldots + 4 f(x_{2N-1}) + f(x_{2N}) \right)\frac {h} {3} $$ত্রুটি#
সিম্পসনের সূত্র দ্বারা সমাকলনের আনুমানিক মানের ত্রুটি হলো
$$ -\tfrac{1}{90} \left(\tfrac{b-a}{2}\right)^5 f^{(4)}(\xi)$$যেখানে $\xi$ হলো $a$ এবং $b$ এর মধ্যে কোনো সংখ্যা।
ত্রুটি অ্যাসিম্পটোটিকভাবে $(b-a)^5$ এর সমানুপাতিক। তবে, উপরের বিশ্লেষণ $(b-a)^4$ এর সমানুপাতিক ত্রুটি নির্দেশ করে। সিম্পসনের রুল একটি অতিরিক্ত অর্ডার লাভ করে কারণ যে বিন্দুগুলোতে ইন্টিগ্র্যান্ড মূল্যায়ন করা হয় সেগুলো $[a, b]$ ব্যবধানে সমানভাবে বিতরিত।
ইমপ্লিমেন্টেশন#
এখানে, $f(x)$ একটি ইউজার-ডিফাইন্ড ফাংশন।
const int N = 1000 * 1000; // number of steps (already multiplied by 2)
double simpson_integration(double a, double b){
double h = (b - a) / N;
double s = f(a) + f(b); // a = x_0 and b = x_2n
for (int i = 1; i <= N - 1; ++i) { // Refer to final Simpson's formula
double x = a + h * i;
s += f(x) * ((i & 1) ? 4 : 2);
}
s *= h / 3;
return s;
}