انتگرال عددی — به زبان ساده (+ دانلود فیلم آموزش رایگان)

در آموزش‌های پیشین مجله فرادرس با مفهوم انتگرال و روش محاسبه انتگرال برخی از توابع مهم، مانند مثلثاتی، هیپربولیک، نمایی، کسری و گنگ آشنا شدیم. برخی توابع هستند که محاسبه مستقیم پادمشتق و انتگرال آن‌ها دشوار است. در چنین مواردی، اگر بخواهیم مقدار یک انتگرال معین را محاسبه کنیم، باید مقدار آن را تقریب بزنیم و انتگرال عددی را به دست آوریم. دو روش بسیار ساده برای این کار وجود دارد که در ادامه آن‌ها را توضیح می‌دهیم. در برخی موارد نیز ناچار به استفاده از روش‌های پیچیده‌تری هستیم.

اگر انتگرال را به عنوان محاسبه مساحت در نظر بگیریم، می‌توان گفت که با جمع مساحت چندین مستطیل مقدار تقریبی انتگرال به دست می‌آید. هرچند این روش برای محاسبه انتگرال عددی ساده است، اما برای آنکه دقت محاسبه مناسب باشد، به تعداد زیادی مستطیل نیاز داریم. یک راه مشابه و البته بهتر، این است که مساحت زیر منحنی را در یک بازه کوچک با مساحت یک ذوزنقه (روش ذوزنقه‌ای) تقریب بزنیم. شکل زیر نشان می‌دهد که با افزایش تعداد زیربازه‌ها دقت تقریب بیشتر شده است.

انتگرال عددی با قاعده ذوزنقه‌ای

برای استفاده از این روش، بازه انتگرال‌گیری را به $$n$$ زیربازه به طول $$\Delta x $$ تقسیم می‌کنیم. یک ذوزنقه معمولی در شکل ۲ نشان داده شده که مساحت آن برابر با $$ { f ( x _ i ) + f ( x _ { i + 1 } ) \over 2 } \Delta x $$ است. اگر مساحت کل ذوزنقه‌ها را با یکدیگر جمع کنیم، خواهیم داشت:

$$ \large \eqalign {
{ f ( x _ 0 ) + f ( x _ 1 ) \over 2 } \Delta x & +{ f ( x _ 1 ) + f ( x_ 2 ) \over 2 } \Delta x + \cdots +
{ f ( x _ { n - 1 } ) + f ( x _ n ) \over2 } \Delta x = \cr
& \left ( { f ( x _ 0 ) \over 2 } + f ( x _ 1 ) + f ( x _2 ) + \cdots + f ( x _ { n - 1 } ) + { f ( x _ n ) \over 2 } \right )
\Delta x . \cr } $$

فیلم آموزش محاسبات عددی – مرور و حل مساله در فرادرس

کلیک کنید

این روش، معمولاً به نام قاعده ذوزنقه‌ای شناخته می‌شود. اگر تعداد زیربازه‌ها خیلی زیاد نباشد، محاسبات با ماشین‌حساب چندان دشوار نخواهد بود. یک کامپیوتر به سادگی برای زیربازه‌های بیشتر قادر به انجام محاسبات است.

در عمل، تقریب انتگرال تنها زمانی مفید است که اندازه دقت را بدانیم. برای مثال، ممکن است به دقتی به اندازه سه رقم اعشار نیاز داشته باشیم. وقتی تقریبی از یک انتگرال را محاسبه می‌کنیم، خطا برابر با اختلاف بین تقریب و مقدار واقعی انتگرال است. از هر روشی که برای انتگرال‌گیری عددی استفاده کنیم، باید یک تخمین خطا (Error Estimate) داشته باشیم. تخمین خطا مقداری است که تضمین شده بزرگ‌تر از خطای واقعی باشد. اگر $$A$$ یک تقریب و $$E$$ تخمین خطای مربوط به آن باشد، آنگاه می‌دانیم که مقدار انتگرال بین $$ A - E $$ و $$ A + E $$ است. در مورد خاص انتگرال‌گیری عددی، می‌خواهیم $$ E = E ( \Delta x ) $$ تابعی برحسب $$ \Delta x $$ باشد که با کوچک شدن $$ \Delta x $$، سریعاً کاهش یابد. خوشبختانه، برای بسیاری از توابع، چنین خطای تخمینی متناظر با تقریب ذوزنقه‌ای وجود دارد.

قضیه: فرض کنید $$ f $$ دارای مشتق دوم $$ f ^ {\prime \prime } $$ در کل بازه $$ [a,b]$$ باشد و برای همه $$x$$ها در این بازه، $$|f^{\prime \prime}(x)|\le M$$. با در نظر گرفتن بازه $$\Delta x= (b-a)/n$$، تخمین خطای تقریب ذوزنقه‌ای به صورت زیر خواهد بود:

$$ \large E ( \Delta x ) = { b - a \over 1 2 } M ( \Delta x ) ^ 2 = { ( b - a ) ^ 3 \over 1 2 n ^ 2 } M . $$

اما چگونه از این قضیه استفاده می‌کنیم؟ مثال زیر پاسخ این پرسش است.

مثال ۱

مقدار انتگرال عددی $$ \int _ 0 ^ 1 e ^ { - x ^ 2 } \, d x $$ را تا دو رقم اعشار بنویسید.

حل: مشتق دوم تابع $$ f=e^{-x^2}$$ برابر با $$ ( 4 x ^ 2 - 2 ) e ^ { - x ^ 2 } $$ است و در بازه $$[0,1]$$، نامساوی $$ | ( 4 x ^ 2 - 2 ) e ^ { - x ^ 2 } | \le 2 $$ را خواهیم داشت. با تخمین تعداد زیربازه‌هایی شروع می‌کنیم که به آن‌ها نیاز داریم. برای داشتن دو رقم اعشار، باید $$E(\Delta x)< 0.005$$ یا

$$ \large \eqalign {
{ 1 \over 1 2 } ( 2 ) { 1 \over n ^ 2 } & < 0 . 0 0 5 \cr
{ 1 \over 6 } ( 2 0 0) & < n ^ 2 \cr
5 . 7 7 \approx \sqrt { 1 0 0 \over 3 } & < n \cr } $$

با $$ n = 6 $$، تخمین خطا $$1/6^3< 0.0047$$ است. اکنون تقریب ذوزنقه‌ای را برای شش زیربازه می‌نویسیم:

$$\large \left ( { f ( 0 ) \over 2 } + f ( 1 / 6 ) + f ( 2 / 6 ) + \cdots + f ( 5 / 6 ) + { f ( 1 ) \over 2 } \right ) { 1 \over 6 }
\approx 0 . 7 4 5 1 2 . $$

بنابراین، مقدار واقعی انتگرال بین $$ 0.74512-0.0047=0.74042$$ و $$ 0.74512+0.0047=0.74982 $$ است. متأسفانه گرد شده مقدار اول $$0.74$$ و گرد شده مقدار دوم $$0.75$$ است. بنابراین، نمی‌توانیم از مقدار صحیح تنها با دو رقم اعشار مطمئن باشیم. در نتیجه، باید یک $$n$$ بزرگ‌تر انتخاب کنیم. همان‌طور که معلوم است، باید $$n=12$$ باشد تا دو محدوده مشابه برای اعشار مورد نظر داشته باشیم. این مقدار $$0.75$$ خواهد بود.

در عمل، معمولاً از مقادیر بهتری برای حداکثر خطای ممکن استفاده می‌شود؛ برای مثال، $$ E ( \Delta x ) < 0 . 0 0 1 $$ انتخاب رایجی است. در این صورت، خواهیم داشت:

$$ \large \eqalign {
{ 1 \over 1 2 } ( 2 ) { 1 \over n ^ 2 } & < 0 . 0 0 1 \cr
{ 1 \over 6 } ( 1 0 0 0 ) & < n ^ 2 \cr
1 2 . 9 1 \approx \sqrt { 5 0 0 \over 3 } & < n \cr } $$

با $$n = 13$$ به نتیجه مطلوب می‌رسیم.

انتگرال عددی با قاعده سیمپسون

تقریب ذوزنقه‌ای عملکرد مناسبی دارد، به ویژه در مقایسه با مستطیل‌ها؛ زیرا وقتی $$ \Delta x $$ به اندازه کافی کوچک باشد، قسمت فوقانی تقریب بهتری از منحنی را ارائه می‌کند. می‌توانیم این ایده را تعمیم دهیم:  اگر بخواهیم تقریبمان دقیق‌تر باشد و به جای یک خط راست از چیز دیگری استفاده کنیم، چه اتفاقی می‌افتد؟ یک گزینه بدیهی برای این کار سهمی است. اگر بتوانیم تکه‌های منحنی را با یک سهمی با معادله $$  y=ax^2+bx+c $$ تقریب بزنیم، می‌توانیم به سادگی سطح زیر نمودار را محاسبه کنیم.

فیلم آموزش محاسبات عددی با متلب MATLAB در فرادرس

کلیک کنید

برای دو نقطه، تعداد سهمی‌هایی که می‌توان تعریف کرد، بینهایت است؛ اما اگر سه نقطه داشته باشیم، تنها یک سهمی وجود دارد که در آن‌ها صدق می‌کند. اگر بخواهیم یک سهمی را پیدا کنیم که از سه نقطه $$ (x_i,f(x_i)) $$، $$(x_{i+1},f(x_{i+1}))$$ و $$ (x_{i+2},f(x_{i+2})) $$ روی منحنی می‌گذرند، در کل بازه $$[x_i,x_{i+2}]$$ بسیار نزدیک به منحنی است. شکل ۳ این موضوع را به خوبی نشان می‌دهد.

اگر بازه $$ [a , b ] $$ را به تعدادی زیربازه زوج تقسیم کنیم، آنگاه می‌توانیم منحنی را با دنباله‌ای از سهمی‌ها تقریب بزنیم که هر کدام دو زیربازه را پوشش می‌دهند. برای اینکه این کار عملی باشد، به فرمول ساده‌ای برای مساحت زیر یک سهمی گذرنده از $$ (x_i,f(x_i))$$، $$(x_{i+1},f(x_{i+1}))$$ و $$(x_{i+2},f(x_{i+2}))$$ نیاز داریم. یعنی باید سهمی را به صورت $$ y=ax^2+bx+c$$ بنویسیم که از این نقاط عبور کند و سپس از آن انتگرال بگیریم. اگرچه جبر مربوط به این کار پیچیده به نظر می‌رسد، اما انجام آن امکان‌پذیر است.

برای یافتن سهمی، سه معادله زیر را برای $$a$$، $$b$$ و $$ c $$ حل می‌کنیم:

$$ \large \eqalign {
f ( x _ i ) & = a ( x _ { i + 1 } - \Delta x ) ^ 2 + b ( x _ { i + 1 } - \Delta x ) + c \cr
f ( x _ { i + 1 } ) & = a ( x _ { i + 1 } ) ^ 2 + b ( x _ { i + 1 } ) + c \cr
f ( x _ { i + 2 } ) & = a ( x _ { i + 1 } + \Delta x ) ^ 2 + b ( x _{ i + 1 } + \Delta x ) + c \cr } $$

پیچیده بودن حل این معادلات غیرمنتظره نیست. با وجود این، می‌توان انتگرال عددی را محاسبه و ساده کرد:

$$ \large \int _ { x _ { i + 1 } - \Delta x } ^ { x _ { i + 1 } + \Delta x } a x ^ 2 + b x + c \, d x =
{ \Delta x \over 3 } ( f ( x _ i ) + 4 f ( x _ { i + 1 } ) + f ( x _ { i + 2 } ) ) . $$

مجموع سطح زیر همه سهمی‌ها برابر است با:

$$ \large \displaylines {
{ \Delta x \over 3 } (f ( x _ 0 ) + 4 f ( x _ { 1 } ) + f ( x _ { 2 } ) + f ( x _ 2 ) + 4 f ( x _ { 3 } ) + f ( x _ { 4 } ) + \cdots
+ f ( x _ { n - 2 } ) + 4 f ( x _ { n - 1 } ) + f ( x _ { n } ) ) = \cr
{ \Delta x \over 3 } ( f ( x _ 0 ) + 4 f ( x _ { 1 }) + 2 f ( x _ { 2 } ) + 4 f ( x _ { 3 } ) + 2 f ( x _ { 4 } ) + \cdots
+ 2 f ( x _ { n - 2 } ) + 4 f ( x _ { n -1 } ) + f ( x _ {n } ) ) . \cr } $$

این فرمول، تنها کمی پیچیده‌تر از فرمول ذوزنقه‌ها است. باید به خاطر داشته باشیم که ضرایب $$2$$ و $$4$$ را قرار دهیم. عدد $$n$$ نیز باید زوج باشد تا معنی داشته باشد. به این تکنیکِ تقریب، قاعده سیمپسون می‌گویند.

مشابه قاعده ذوزنقه‌ای، یک تخمین خطا برای قاعده سیمپسون نیز وجود دارد.

قضیه: فرض کنید $$f$$ دارای مشتق چهارم $$ f ^ {(4)}$$ در کل بازه $$[a,b]$$ بوده و برای همه $$x$$های این بازه، $$|f^{(4)}(x)|\le M$$. اگر $$\Delta x= (b-a)/n$$، یک تخمین خطا برای تقریب سیمپسون به صورت زیر خواهد بود:‌

$$ \large E ( \Delta x ) = { b - a \over 1 8 0 } M ( \Delta x ) ^ 4 ={ ( b - a ) ^ 5 \over 1 8 0 n ^ 4 } M . $$

مثال زیر، کاربرد این قضیه را نشان می‌دهد.

مثال ۲

انتگرال عددی $$ \int_0^1 e^{-x^2}\,dx $$ را تا دو رقم اعشار تقریب بزنید.

حل: مشتق چهارم $$ f=e^{-x^2}$$ برابر با $$(16x^2-48x^2+12)e^{-x^2}$$ بوده و روی بازه $$ [0,1]$$ حداکثر مقدار مطلق آن برابر با $$12$$ است. با تخمین تعداد زیربازه‌هایی که نیاز داریم شروع می‌کنیم. برای داشتن دقت دو رقم اعشار، باید $$ E(\Delta x)< 0.005 $$ را داشته باشیم. اما به توجه به مثال قبل، لازم است $$ E(\Delta x)< 0.001 $$ را داشته باشیم:

بنابرین، $$ n = 4 $$ را امتحان می‌کنیم، زیرا به تعداد زوجی زیربازه نیاز داریم. سپس تخمین خطا $$ 12/180/4^4< 0.0003 $$ خواهد بود و تقریب برابر است با:

$$ \large ( f ( 0 ) + 4 f ( 1 / 4 ) + 2 f ( 1 / 2) + 4 f ( 3 / 4 ) + f (1 ) ) { 1 \over 3 \cdot 4 }
\approx 0 . 7 4 6 8 5 5 . $$

بنابراین، مقدار واقعی انتگرال بین $$ 0.746855-0.0003=0.746555 $$ و $$ 0.746855+0.0003=0.7471555 $$ خواهد بود که هر دوی آن‌ها به $$ 0.75$$ گرد می‌شوند.

انتگرال عددی با قاعده سیمپسون در متلب

برنامه متلب قاعده سیمپسون نیز به صورت زیر است:

فیلم آموزش محاسبات عددی با متلب MATLAB در فرادرس

کلیک کنید

1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2% Simpson Rule (Version 1.0)                            %
3%                                                                %
4% Programmed By: MatlabSite.com Team.                            %
5%                Copyright 2009                                  %
6%                                                                %
7%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8! Copyright.pdf
9
10
11clc;
12clear;
13x=sym('x');
14
15disp('This program calculates integral of f(x) using Composite Simpson rule.');
16f=input('Enter f(x) as a function of x: ');
17disp(' ');
18
19a=input('Enter the start of integration range: ');
20b=input('Enter the end of integration range: ');
21disp(' ');
22
23xmin=min(a,b);
24xmax=max(a,b);
25
26a=xmin;
27b=xmax;
28
29n=input('Enter the number of subintervals (an even number): ');
30disp(' ');
31n=2*floor(n/2);
32if(n<2)
33    n=2;
34end
35
36h=(b-a)/n;
37S=0;
38for i=1:n+1
39    xx=a+(i-1)*h;
40    ff=subs(f,xx);
41    
42    if i==1 || i==n+1
43        S=S+ff;
44    else
45        if mod(i,2)==0
46            S=S+4*ff;
47        else
48            S=S+2*ff;
49        end
50    end
51end
52
53S=S*h/3;
54
55disp(['Integral of ' char(f) ' in the range [' mat2str(a,4) ', ' mat2str(b,4) '] is:']);
56disp(mat2str(S,10));
57disp(' ');

انتگرال عددی با قاعده ذوزنقه‌ای در متلب

برنامه متلب قاعده ذوزنقه‌ای نیز به صورت زیر است:

فیلم آموزش محاسبات عددی با متلب MATLAB در فرادرس

کلیک کنید

1function I = trapezoidal_f1 ( f )
2% for calculating integrals using trapezoidal rule when function is known
3
4%asking for the range and desired accuracy
5 R= input( ' Enter the limits of integrations [ x_min, x_max]   :\n');
6 tol = input(' Error allowed in the final answer should be of an order : \n');
7 a= R(1,1); b = R(1,2);
8 
9 %initial h and n
10 n = 100;
11 h = (b -a )/100;
12 
13 %for calculating maximum of h^2 *f''(x) in the given region
14 for k = 0:100
15        x( 1, k+1 ) = a +  k*h ;
16        y2 ( 1, k+1) = feval ( f, x(1,k+1) + 2*h ) - 2*feval( f, x(1,k+1) + h ) + feval( f, x(k+1) );
17 end 
18 [ y i ] =  max( y2);
19 x_opt = x(1,i);
20   
21 % for calculating the desired value of h
22 m=0;
23 while abs((feval ( f, x_opt + 2*h ) - 2*feval(f,  x_opt+ h) + feval(f, x_opt)) * ( b-a)/12) > tol % Global error for trapezoidal rule
24                m = m +1;
25                h = h * 10^-m;
26 end
27 
28 %calculating  n
29 n = ceil( (b-a)/h );
30 h = ( b-a )/ n;
31
32 % forming matrix X
33 for k = 1:(n+1)
34    X(k,1) = a + (k-1)*h;
35    X(k,2) = feval( f, X(k,1));
36 end
37 
38% trapezoidal formula
39I = h/2 * (2*sum( X(:,2))- X(1,2)- X( n,2));
40 % displaying final result
41 h 
42 n
43 disp(sprintf(' The area under this curve in the interval ( %10f , %10f )  is %10.6f .',a,b,I))

اگر این مطلب برایتان مفید بوده است، آموزش‌های زیر نیز به شما پیشنهاد می‌شوند:

• مجموعه آموزش‌های دروس ریاضیات
• آموزش محاسبات عددی با MATLAB
• مجموعه آموزش‌های محاسبات عددی
• آموزش ریاضی پایه دانشگاهی
• روش نیوتن — به زبان ساده
• دستگاه معادلات خطی — به زبان ساده
• جذر یا محاسبه ریشه دوم عدد — به زبان ساده

^^