پرش به محتویات

انتگرال‌گیری به روش سیمپسون

می‌خواهیم مقدار یک انتگرال معین را محاسبه کنیم

$$\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)$ را با یک سهمی $P(x)$ که از ۳ نقطه $(x_ {2i-2}, x_ {2i-1}, x_ {2i})$ می‌گذرد، جایگزین می‌کنیم. چنین سهمی همیشه وجود دارد و یکتا است؛ می‌توان آن را به روش تحلیلی پیدا کرد. برای مثال، می‌توانیم آن را با استفاده از درون‌یابی چندجمله‌ای لاگرانژ (Lagrange polynomial interpolation) بسازیم. تنها کاری که باقی می‌ماند، انتگرال‌گیری از این چندجمله‌ای است. اگر این کار را برای یک تابع عمومی $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; // تعداد گام‌ها (قبلاً در ۲ ضرب شده است)

double simpson_integration(double a, double b){
    double h = (b - a) / N;
    double s = f(a) + f(b); // a = x_0 و b = x_2n
    for (int i = 1; i <= N - 1; ++i) { // به فرمول نهایی سیمپسون مراجعه کنید
        double x = a + h * i;
        s += f(x) * ((i & 1) ? 4 : 2);
    }
    s *= h / 3;
    return s;
}

مسائل تمرینی