برون یابی ریچاردسون — از صفر تا صد

در این آموزش از مجموعه مطالب ریاضی مجله فرادرس، با «برون یابی ریچاردسون» (Richardson Extrapolation) آشنا می‌شویم و مشتق‌گیری و انتگرال‌گیری با استفاده از این روش را توضیح خواهیم داد.

برون یابی ریچاردسون

دقت محدوده وسیعی از روش‌های عددی برای تقریب تکراری مقدار متغیر مشخص $$A$$ به اندازه گام $$h$$ وابسته است. در حالت کلی، کاهش $$h$$ مرتبه خطای برش $$O ( h ^ k )$$ نتیجه را افزایش خواهد داد. از لحاظ نظری، مقدار دقیق $$A$$ را می‌توان با حد $$h \to \infty$$ به دست آورد. برون یابی ریچاردسون یک روش عمومی است که عملکرد تکرار را با کاهش اندازه گام و در نتیجه افزایش مرتبه خطای «برش» (Truncation) را بهبود می‌دهد.

فیلم آموزش آنالیز عددی پیشرفته – جامع و با نکات مهم در فرادرس

کلیک کنید

نماد $$A _ 0 ( h )$$ را به عنوان تقریب $$A$$ در نظر می‌گیریم که با یک روش مشخص بر اساس اندازه گام $$h$$ به دست می‌آید و فرض می‌کنیم مرتبه خطای برش $$O ( h ^ {k_1} )$$ باشد:

$$\large \begin {equation} A = A _ 0 ( h ) + C _ 1 h ^ { k _1 } + C _ 2 h ^ { k _ 2 } + C _ 3 h ^ { k _ 3 } + \cdots = A _ 0 ( h ) + O ( h ^ { k _ 1 } ) \end {equation} \quad \quad (1)$$

که $$k_1$$. با جایگزینی اندازه گام $$h$$ با $$h / t$$ برای $$t > 1$$، خواهیم داشت:

$$\large \begin {equation} A = A _ 0 ( h / t ) + C _ 1 ( h / t ) ^ { k _ 1 } + C _ 2 ( h / t ) ^ { k _ 2 } + \cdots = A _ 0 ( h / t ) + O ( h / t ) ^ { k _ 1 } \end {equation} \quad \quad (2)$$

با کم کردن معادله اول از $$t ^ { k_ 1 }$$ برابر معادله دوم، داریم:

$$\large {t ^ { k _ 1 } A - A} = t ^ { k _ 1 } \left [ A _ 0 ( h / t ) + C _ 1 ( h / t ) ^ { k _ 1 } + C _ 2 ( h / t ) ^ { k _ 2 } + \cdots ] \\ - [ A _ 0 ( h ) + C _ 1 h ^ { k _ 1 } + C _ 2 h ^ { k _ 2 } + \cdots \right ] \quad \quad ( 3 )$$

و در نتیجه، خواهیم داشت:‌

$$\large \begin {equation} A = \frac { t ^ { k _ 1 } A _ 0 ( h / t ) - A _ 0 ( h ) }{ t ^ { k _ 1 } - 1 } + O ( h ^ { k _ 2 } ) = A _ 1 ( h / t ) + O ( h ^ { k _ 2 } ) \end {equation} \quad \quad ( 4 )$$

بنابراین، یک روش جدید با گام زمانی $$h / t$$ داریم:‌

$$\large \begin {equation} A _ 1 ( h / t ) = \frac { t ^ { k _ 1 } A _ 0 ( h / t ) - A _ 0 ( h ) } { t ^ { k _ 1 } - 1 } \end {equation} \quad \quad ( 5 )$$

با یک جمله خطای کوچک‌تر از مرتبه $$O ( h ^ { k _ 2 } )$$، بزرگ‌تر از $$O ( h ^ {k_1})$$ روش قدیمی $$A _ 0 ( h )$$. این روش برون یابی ریچاردسون را می‌توان برای به دست آوردن روش $$A _ 2$$ دیگر با یک جمله خطای مرتبه بالاتر $$O(h^{k^3})$$، به صورت بازگشتی به $$A _ 1$$ اعمال کرد.

بازگشت برون یابی ریچاردسون به صورت زیر نشان داده شده است. روش $$A _ 0 ( h )$$ را که مقدار $$A$$ را بر اساس اندازه گام $$h$$ تقریب می‌زند، در نظر بگیرید:

$$\large \begin {equation} A = A _ 0 ( h ) + C _ 1 h ^ 2 + C _ 2 h ^ 4 + C _ 3 h ^ 6 + C _ 4 h ^ 8 + \cdots = A _ 0 ( h ) + O ( h ^ 2 ) \end {equation} \quad \quad (6 )$$

اکنون از برون یابی ریچاردسون با $$t = 2$$ برای بهبود دقت تقریب استفاده می‌کنیم:

• $$n = 1$$: با جایگزینی $$h$$ در معادله با $$h / 2$$، داریم:

$$\large \begin {equation} A = A _ 0 ( h / 2 ) + C _ 1 \frac { h ^ 2 } { 2 ^ 2} + C _ 2 \frac { h ^ 4 } { 2 ^ 4 } + C _ 3 \frac { h ^ 6 } { 2 ^ 6 } + C _ 4 \frac { h ^ 8 } { 2 ^ 8 } \cdots , \; \; \; \end {equation} \quad \quad (7)$$

با ضرب دو طرف در $$t ^ k = 2 ^ 2 = 4$$، خواهیم داشت:

$$\large \begin {equation} 4 A = 4 A _ 0 ( h / 2 ) + C _ 1 h ^ 2 + C _ 2 \frac { h ^ 4 }{ 4 } + C _ 3 \frac { h ^ 6 } { 1 6 } + C _ 4 \frac { h ^ 8 } { 6 4 } \cdots \end {equation} \quad ( 8 )$$

و با کم کردن معادله (۶)، خواهیم داشت:

$$ \large \begin {eqnarray}<br /> A & = & \frac { 4 A _ 0 ( h / 2 ) - A _ 0 ( h ) } { 3 } -C _ 2 \frac { h ^ 4 } { 4 } - C _ 3 \frac { 5 h ^ 6 } { 1 6 }<br /> - C _ 4 \frac { 2 1 h ^ 8 } { 6 4 } + \cdots<br /> \nonumber \\<br /> & = & A _ 1 ( h / 2 ) - C _ 2 \frac { h ^ 4 } { 4 } -C _ 3 \frac { 5 h ^ 6 } { 1 6 } - C _ 4 \frac { 2 1 h ^ 8 } { 6 4 } + \cdots<br /> = A _ 1 ( h / 2 ) + O ( h ^ 4 )<br /> \end {eqnarray} \quad (9) $$

در اینجا، یک روش جدید را تعریف می‌کنیم:

$$\large \begin {equation} A _ 1 ( h / 2 ) = \frac { 4 A _ 0 ( h / 2 ) - A _ 0 ( h ) } { 3 } \end {equation} \;\;\;\; ( 10 )$$

با یک اندازه گام $$h / 2$$ و جمله خطای $$O ( h ^ 4 )$$.

• $$n = 2$$: با جایگزینی $$h$$ در معادله (۹) با $$h / 2$$، مجدداً خواهیم داشت:

$$\large \begin {equation} A = A _ 1 ( h / 4 ) - C _ 2 \frac { h ^ 4 } { 4 \times 2 ^ 4 } -C _ 3 \frac { 5 h ^ 6 } { 1 6 \times 2 ^ 6 } - C _ 4 \frac { 2 1 h ^ 8 } { 6 4 \times 2 ^ 8 } + \cdots \end {equation} \; \;\;\; ( 11)$$

با ضرب دو طرف در $$2 ^ 4 = 16$$، می‌توان نوشت:

$$\large \begin {equation} 16 A = 1 6 A _ 1 ( h / 4 ) - C _ 2 \frac { h ^ 4 } { 4 } -C _ 3 \frac { 5 h ^ 6 } { 6 4 } - C _ 4 \frac { 2 1 h ^ 8 } { 1 0 2 4 } + \cdots \end {equation} \;\;\; (12)$$

و با کم کردن معادله (۹)، داریم:

$$ \large \begin {eqnarray}<br /> A & = & \frac { 1 6 A _ 1 ( h / 4 ) - A _ 1 ( h / 2 ) }{ 1 5 } + C _ 3 \frac { h ^ 6 } { 6 4 }<br /> + C _ 4 \frac { 2 1 h ^ 8 } { 1 0 2 4 } + \cdots<br /> \nonumber \\<br /> & = & A _ 2 ( h / 4 ) + C _ 3 \frac { h ^ 6 }{ 6 4 } + C _ 4 \frac { 2 1 h ^ 8 } { 1 0 2 4 } + \cdots<br /> = A _ 2 ( h / 4 ) + O ( h ^ 6 )<br /> \end {eqnarray} \;\;\; \; ( 1 3 ) $$

در اینجا، یک روش جدید تعریف شده است:

$$\large \begin {equation} A _ 2 ( h / 4 ) = \frac { 1 6 A _ 1 ( h / 4 ) - A _ 1 ( h / 2 ) } { 1 5 } \end {equation} \;\;\; (14)$$

• $$n=3$$: با جایگزینی $$h$$ در معادله (۱۱)، با $$h / 2$$، مجدداً خواهیم داشت:

$$\large \begin {equation} A = A _ 2 ( h / 8 ) + C _ 3 \frac { h ^ 6 } { 6 4 \times 2 ^ 6 } + C _ 4 \frac { 2 1 h ^ 8 } { 1 0 2 4 \times 2 ^ 8 } + \cdots \end {equation} \;\;\; ( 15 )$$

با ضرب دو طرف در $$2 ^ 6 = 64$$، خواهیم داشت:

$$\large \begin {equation} 6 4 A = 6 4 A _ 2 ( h / 8 ) + C _ 3 \frac { h ^ 6 }{ 6 4 } + C _ 4 \frac { 2 1 h ^ 8 } { 4 0 9 6 } + \cdots \end {equation} \;\;\; (16)$$

و با کم کردن معادله (۱۳)، می‌توان نوشت:

$$ \large \begin {eqnarray}<br /> A & = & \frac { 6 4 A _ 2 ( h / 8 ) - A _ 2 ( h / 4 ) } { 6 3 } - C _ 4 \frac { h ^ 8 } { 4 0 9 6 } + \cdots<br /> \nonumber \\<br /> & = & A _ 3 ( h / 8 ) + O ( h ^ 8 )<br /> \end {eqnarray} \;\;\; (17)$$

در اینجا یک روش جدید دیگر داریم:

$$\large \begin {equation} A _ 3 ( h / 8 ) = \frac { 6 4 A _ 2 ( h / 8 ) - A _ 2 ( h / 4 ) }{ 6 3 } \end {equation} \;\;\;\; ( 18 )$$

با یک اندازه گام $$h / 8$$ و جمله خطای $$O ( h ^ 8 )$$.

می‌بینیم که برای هر تکرار در فرایند، مرتبه جمله خطا با ۲ افزایش می‌یابد. اکنون به جای $$A _ n ( h / 2 ^ m )$$ از نماد $$A _ {m,n}$$ استفاده می‌کنیم که $$m$$ برای اندازه گام $$h _ m = h _ { m -1} / 2 = \cdots  = h _ 0 / 2 ^ m$$ با $$h _ 0 = h$$ است، و $$n$$ برای روش $$A _ n$$ در $$n$$اُمین سطح بازگشت ($$m \ge n$$) به دست آمده و فرمول کلی بازگشتی زیر را نشان می‌دهد:

$$\large \begin {equation} A _ { m , n } = \frac { 4 ^ n A _ { m , n - 1 } - A _ { m - 1 , n - 1 } } { 4 ^ n - 1 } \end {equation} (19)$$

این فرایند بازگشتی را می‌توان به فرم جدولی زیر نوشت:

در اینجا، $$m$$ اندیس بعد عمودی است که اندازه‌های گام $$h _ m = h / 2 ^ m$$ را نشان می‌دهد، و $$n$$ اندیس بعد افقی است که روش به دست آمده در سطح $$n$$ بازگشت را نشان می‌دهد. عنصر عمومی $$A _ { m , n }$$ را می‌توان با رابطه بازگشتی در معادله (۱۹) بر اساس همسایه سمت چپ آن $$A _{ m , n- 1 }$$ بر اساس روش قبلی $$A _ { n - 1 }$$ با اندازه گام مشابه $$h _ m$$ به دست آورد. همسایه بالا چپ $$A _ { m - 1 , n - 1 }$$ نیز با روش قبلی $$A _ { n - 1 }$$ تولید شده است، اما با استفاده از یک اندازه گام $$h _ { m - 1 } = 2 h _ m$$. جمله خطای $$A _ { n , n }$$ روی قطر جدول دو درجه بزرگ‌تر از $$A _ { n - 1 , n - 1 }$$ و $$2 n$$ درجه بزرگ‌تر از روش اولیه $$A _ { 0 , 0 }$$ است.

مشتق گیری با برون یابی ریچاردسون

برون یابی ریچاردسون را می‌توان برای بهبود مشتق‌گیری عددی به کار برد. مجدداً ذکر می‌کنیم که $$D [ f ] = f' ( x )$$ را می‌توان با تفاضل مرکزی مبتنی بر دو نقطه همسایه $$x - h$$ و $$x + h$$ با اندازه گام $$h$$ و جمله خطای برش $$O ( h ^ 2 )$$، مطابق معادله زیر به دست آورد:

$$ \large \begin {eqnarray}<br /> D [ f ] & = & \frac { f ( x + h ) - f ( x - h ) } { 2 h }<br /> - \left [ \frac { f ^ {\prime \prime \prime } ( x ) } { 3 ! } h ^ 2 + \frac { f ^ { ( 5 ) } ( x ) } { 5 ! } h ^ 4 + \cdots \right ]<br /> \nonumber \\ & = & D _ 0 ( h ) + C _ 1 h ^ 2 + C _ 2h ^ 4 + C _ 3 h ^ 6 + C _ 4 h ^ 8 + \cdots = D _ 0 ( h ) +C _2 h ^ 2 + O ( h ^ 4 )<br /> \end {eqnarray} \;\;\; (20)$$

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

کلیک کنید

که در آن:

$$\large \begin {equation} D _ 0 ( h ) = \frac { f ( x + h ) - f ( x - h ) } { 2 h } \end {equation} \;\;\; (21)$$

دقت این روش $$D _ 0 ( h )$$ را می‌توان از قبل با برون یابی ریچاردسون بهبود داد. ابتدا $$h$$ را در معادله با $$h / 2$$ جایگزین می‌کنیم:

$$\large \begin {equation} D [ f ] = D _ 0 ( h / 2 ) + C _ 2 h ^ 2 / 4 + O ( h ^ 4 ) \end {equation} \;\;\; (22)$$

و آنگاه، خواهیم داشت:

$$\large \begin {equation} D [ f ] = \frac { 4 D _ 0 ( h/ 2 ) - D _ 0 ( h ) }{ 3 } + O ( h ^ 4 ) = D _ 1 ( h / 2 ) + O ( h ^ 4 ) \end {equation} \;\;\; (23)$$

که در آن،

$$ \large \begin {eqnarray}<br /> D _ 1 ( h / 2 ) & = & \frac { 4 D _ 0 ( h / 2 ) - D _ 0 ( h ) } { 3 }<br /> = \frac { 1 } { 3 } \left [ 4 \frac { f ( x + h / 2 ) - f ( x - h / 2 ) } { h }<br /> - \frac { f ( x + h ) - f ( x - h ) } {2 h } \right ]<br /> \nonumber \\<br /> & = & \frac { 1 } { 6 h } \left [ f ( x - h ) - 8 f ( x - h / 2 ) + 8 f ( x + h / 2 ) - f ( x + h ) \right ]<br /> \end {eqnarray} \;\;\; ( 24) $$

یک روش جدید مبتنی بر چهار نقطه $$x - h$$، $$x - h / 2$$، $$x + h / 2$$ و $$x + h$$، با اندازه گام $$h / 2$$ و یک جمله خطای $$O ( h ^ 4 )$$ است.

دقت را می‌توان با فرایند مشابهی بیشتر بهبود داد. مجدداً $$h$$ را در $$D _ 1 ( h / 2 )$$ با $$h / 2$$ جایگزین کرده و گام‌های مشابهی را طی می‌کنیم:

$$\large \begin {equation} D [ f ] = \frac { 1 6 D _ 1 ( h / 4 ) - D _ 1 ( h / 2 ) }{ 1 5 } + O ( h ^ 6 ) = D _ 2 ( h / 4 ) + O ( h ^ 6 ) \end {equation} \;\;\; ( 25 )$$

که بر اساس شش نقطه $$x \pm h$$، $$x \pm h/ 2$$ و $$x \pm h/ 4$$ با اندازه گام $$h / 4$$ و جمله خطای $$O ( h ^ 6 )$$ منجر به یک روش جدید می‌شود:

$$ \large \begin {eqnarray}<br /> D _ 2 ( h /4 ) & = & \frac { 1 6 D _ 1 ( h / 4 ) - D _ 1 ( h / 2 ) } { 1 5 }<br /> \nonumber \\<br /> & = & \frac { 1 } { 1 5 } \; \frac { 1 6 } { 3 h } [ f ( x -h / 2 ) - 8 f ( x - h / 4 )<br /> + 8 f ( x + h / 4 ) + f ( x + h / 2 ) ]<br /> \nonumber \\<br /> & & - \frac { 1 } { 1 5 } \; \frac { 1 } { 6 h } [ f ( x - h ) - 8 f ( x - h / 2 ) + 8 f ( x + h / 2 ) - f ( x + h ) ]<br /> \nonumber \\<br /> & = & \frac { 1 } { 9 0 h } [ - f ( x - h ) + 4 0 f ( x -h / 2 ) - 2 5 6 f ( x - h / 4 )<br /> \nonumber \\<br /> & & + 2 5 6 f ( x + h / 4 ) - 4 0 f ( x + h / 2 ) + f ( x + h ) ]<br /> \end {eqnarray} \;\;\; (26) $$

لازم به ذکر است که $$D _ 1 ( h )$$ و $$D _ 2 ( h )$$ در واقع روش‌های مشابه قبل هستند که با روش ضرایب نامعین به دست آمدند. با این حال، اکنون می‌دانیم فرایندی را که در بالا مورد بحث قرار گرفت، می‌توان به صورت بازگشتی و برای به دست آوردن روش‌های $$D_ 3$$، $$D_ 4$$ و... با نتایج دقیق‌تر انجام داد.

فرایند برون یابی ریچاردسون را می‌توان به صورت تکراری برای بهبود دقت انجام داد.

مثال مشتق گیری با برون یابی ریچاردسون

می‌خواهیم مشتق تابع گاوسی $$f ( x ) = e ^ { - x ^ 2 }$$ در $$x = 1$$ را به دست آوریم. ابتدا، $$n = m = 0$$ را قرار می‌دهیم و از روش تفاضل تقسیم شده با اندازه گام $$h _ 0 = h = 1$$ استفاده می‌کنیم:‌

$$\large \begin {equation} f' ( x ) \approx D _ { 0 , 0 } = \frac { f ( x + h ) - f ( x - h ) }{ 2 h } = - 0 . 4 9 0 8 \end {equation} \;\;\; (27)$$

 از $$n = 0$$ برای روش مشابه شروع می‌کنیم، اما $$m = 1$$ را برای یک گام کوچک‌تر سمت $$h _ 1 = h _ 0 / 2 ^ m  = h _ 0 / 2 = 0.5$$:

$$\large \begin {equation} D _ { 1 , 0 } = \frac { f ( x + h / 2 ) - f ( x - h / 2 ) }{ 2 } = - 0 . 6 7 3 4 \end {equation} \;\;\; ( 28)$$

اکنون یک روش جدید با استفاده از برون یابی ریچاردسون برای ایجاد یک نتیجه دقیق‌تر به دست می‌آوریم:

$$\large \begin {equation} D _ { 1 , 1 } = \frac { 4 D _ { 1 , 0 } - D _ { 0 , 0 } }{ 3 } = - 0 . 7 3 4 2 5 \end {equation} \;\; (29)$$

این فرایند به صورت بازگشتی اغلب برای تولید نتایج نشان داده شده در جدول زیر انجام می‌شود. هر وقت اندیس سطر $$m$$ به اندازه ۱ اضافه شود، اندازه گام $$h_m=h/2^m$$ نصف می‌شود؛ و هر زمان اندیس ستون $$n$$ به اندازه  زیاد شود، یک روش جدید با یک مرحله بیشتر از بازگشت به دست می‌آید. همه نتایج در $$n$$اُمین ستون با روش مشابه $$D _ n$$ اما بر اساس اندازه گام‌های مختلف تولید می‌شود، اما همه بر اساس اندازه گام مشابه $$h_n=h/2^n$$. هر عنصر $$D_{m,n}$$ در جدول بر اساس همسایه چپ $$D_{m,n-1}$$ و همسایه بالا چپ $$D_{m-1,n-1}$$ با بازگشت به دست می‌آید:

$$\large \begin {equation} D _ { m , n } = \frac { 4 ^ n D _ { m , n - 1 } - D _ { m - 1 , n - 1 } } { 4 ^ n - 1 } \end {equation} \;\; (30)$$

دقیق‌ترین نتیجه در هر مرحله بازگشت روی قطر جدول است. در این مورد خاص، $$D_{4,4}=-0.7357589$$ دقیق‌ترین نتیجه است که بعد از $$n = 4$$ سطح بازگشت با استفاده از یک اندازه گام $$h_4=h/2^4=h/16$$ به دست می‌آید.

سطر آخر نشان می‌دهد که خطاهای نسبی پنج روش $$D_{n,n}, \;n=0,\cdots,4$$ در قطر اصلی با مقدار واقعی مقایسه شده‌اند:

$$\large \begin {equation} f' ( 1 ) = \frac { d } { d x } e ^ { - x ^ 2 } \bigg | _ { x = 1 } = - 2 x \, e ^ { - x ^ 2 } \bigg | _ { x = 1 } = - 0.73575888 \end {equation} \;\;\; (31)$$

کد متلب زیر این جدول را نتیجه داده است:

1    h=1;                                    % inital step size
2    n=5;                                    % number of recursions
3    m=n;
4    D=zeros(n);                             % initializing array D
5        for i=1:m                           % for all rows 
6            n=2^(i-1);
7            D(i,1)=(f(x+h)-f(x-h))/2/h;
8            for j=2:i                       % for all columns
9                D(i,j)=(4^(j-1)*D(i,j-1)-D(i-1,j-1))/(4^(j-1)-1);
10            end
11            h=h/2;
12        end

انتگرال گیری با برون یابی ریچاردسون

برون یابی ریچاردسون را می‌توان برای بهبود انتگرال‌گیری عددی نیز به کار برد (با روش منتجه‌ای که «روش رمبرگ» (Romberg Method) نامیده می‌شود). با یادآوری قانون ذوزنقه‌ای، خواهیم داشت:

$$ \large \begin {eqnarray}<br /> I [ f ] & = & \int _ a ^ b f ( x ) \, d x<br /> \approx h \frac { f ( a ) + f ( b ) } { 2 } + h \, \sum _ { k = 1 } ^ { n - 1 } f ( a + k h )<br /> \nonumber \\<br /> & = & h \frac { f ( a ) - f ( b ) } { 2 } + h \,\sum _ { k = 1 } ^ n f ( a + k h )<br /> \end {eqnarray} \;\;\; (32) $$

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

کلیک کنید

که $$h = ( b - a ) / n$$ است. در تئوری، با حد $$n \to \infty$$ خواهیم داشت: $$( b - a ) / n \to 0$$ و $$I _ 0 ( h )$$ انتگرال $$I [ f ]$$ خواهد بود:

$$\large \begin {equation} \lim \limits _ { h \rightarrow 0 } I _ 0 ( h ) = \int _ a ^ b f ( x ) \, d x = I [ f ] \end {equation} \;\; (33)$$

البته این روش از نظر محاسباتی عملی نیست.

با این حال، می‌توانیم از برون یابی ریچاردسون استفاده کنیم تا درجه جملات خطا را افزایش داده و دقت را بهبود بخشیم. به طور خاص، یک متغیر دیگر $$t$$ را معرفی می‌کنیم، به گونه‌ای که $$x=a+t\,h$$، $$dx=h\,dt$$، $$x _ k = a + kh$$ و $$b = a + n h$$. علاوه بر این، $$\phi(t)=f(a+t\,h)$$ را تعریف می‌کنیم و خواهیم داشت:

$$\large \begin {equation} \phi ^ { ( k ) } ( t ) = \frac { d ^ k } { d t ^ k } \phi ( t ) = \frac { d ^ k } { d t ^ k } f ( a + t h ) = h ^ k f ^ { ( k ) } ( x ) \end {equation} \;\;\; (34)$$

اکنون $$f ( a )$$، $$f ( a + kh)$$ و $$f ( b ) = f ( a + n h)$$ به ترتیب، به $$\phi ( 0 )$$، $$\phi (k)$$ و $$\phi (n )$$ تبدیل خواهند شد و روش ذوزنقه‌ای بالا به صورت زیر در می‌آید:

$$ \large \begin {eqnarray}<br /> I [ f ] & = & \int _ a ^ b f ( x ) \, d x = h \int _ 0 ^ n \phi ( t ) \, d t \approx<br /> h \frac { \phi ( 0 ) - \phi ( n ) } { 2 } + h \sum _ { k = 1 } ^ n \phi ( k )<br /> \end {eqnarray} \;\;\; (35) $$

با مقایسه این با فرمول اویلر-مکلورن، داریم:

$$\large \begin {equation} \int _ 0 ^ n \phi ( t ) \, d t = \sum _ { k = 1 } ^ n \phi ( k ) + \frac { \phi ( 0 ) - \phi ( n ) } { 2 } - \sum _ { k = 1 } ^ \infty \frac { B _ { 2 k } } { ( 2 k ) ! }[ \phi ^ { ( 2 k - 1 ) } ( n ) - \phi ^ { ( 2 k - 1 ) } ( 0 ) ] \end {equation} \;\;\; (36)$$

درنتیجه، خواهیم داشت:

$$ \large \begin {eqnarray}<br /> I [ f ] & = & \int _ a ^ b f ( x ) \, d x = h \int _ 0 ^ n \phi ( t ) \, d t<br /> \nonumber \\<br /> & = & h \frac { \phi ( 0 ) - \phi ( n ) }{ 2 } + h \sum _ { k = 1 } ^ n \phi ( k )<br /> - h \sum _ { k = 1 } ^ \infty \frac { B _ { 2 k } } { ( 2 k ) ! } [ \phi ^ { ( 2 k - 1 ) } ( n ) - \phi ^ { ( 2 k - 1 ) } ( 0 ) ]<br /> \nonumber \\<br /> & = & h \left [ \frac { f ( a ) - f ( b ) }{ 2 } + \sum _ { k = 1 } ^ n f ( a + k h ) \right ]<br /> - h ^ { 2 k } \sum _ { k = 1 } ^ \infty \frac { B _ { 2 k } }{ ( 2 k ) ! } [ f ^ { ( 2 k - 1 ) } ( b ) - f ^ { ( 2 k - 1 ) } ( a ) ]<br /> \nonumber \\<br /> & = & I _ 0 ( h ) + \sum _ { k = 1 } ^ \infty C _ k h ^ { 2 k } = I _ 0 ( h ) + C _ 1 h ^ 2 + C _ 4 h ^ 4 + \cdots<br /> = I _ 0 ( h ) + O ( h ^ 2 )<br /> \end {eqnarray} \;\;\; (37) $$

که $$I _ 0 ( h )$$ روش قاعده ذوزنقه‌ای ترکیبی را مشخص می‌کند:

$$\large \begin {equation} I _ 0 ( h ) = h \frac { f ( a ) - f ( b ) } { 2 } + h \sum _ { k = 1 } ^ n f ( a + k h ) \end {equation} \;\;\; ( 38 )$$

با اندازه گام $$h = (b - a )/ n$$ و جمله خطای $$O ( h ^ 2 )$$، و

$$\large \begin {equation} C _ k = - \frac { B _ { 2 k } } { ( 2 k ) ! } [ f ^ { ( 2 k - 1 ) }( b ) - f ^ { ( 2 k - 1) } ( a ) ] , \; \; \; \; \; ( k = 1 , \, 2 , \, 3 , \cdots ) \end {equation} \;\;\; (39)$$

اکنون می‌توانیم از برون یابی ریچاردسون برای بهبود دقت این روش $$I _ 0 ( h )$$، با تعویض $$h$$ معادله (37) با $$h / 2$$ استفاده کنیم:

$$\large \begin {equation} I [ f ] = I ( h / 2 ) + C _ 1 \left ( \frac { h } { 2 } \right ) ^ 2 + C _ 2 \left ( \frac { h } { 2 } \right ) ^ 4 + O ( h ^ 6 ) \end {equation} \;\;\; (40)$$

و آنگاه معادله زیر را از ۴ برابر این معادله کم می‌کنیم:

$$\large \begin{equation} \left[\begin{array}{rcc}1 & 1 & 1\\-1 & 0 & 1\\1 & 0 & 1\end{array}\right] [{\bf c}_1,\,{\bf c}_2] =\left[\begin{array}{cc}0&0\\1/h&0\\0&2/h^2\end{array}\right] \end{equation} \;\;\; (41)$$

و خواهیم داشت:

$$\large \begin {equation} I [ f ] = \frac { 4 I ( h / 2 ) - I _ n ( h ) } { 3 } + O ( h ^ 4 ) = I _ 1 ( h / 2 ) + O ( h ^ 4 ) \end {equation} \;\;\; (42)$$

که در آن،

$$\large \begin {equation} I _ 1 ( h / 2 ) = \frac { 4 I ( h / 2 ) - I ( h ) } { 3 } \end {equation} \;\;\;\; ( 43 )$$

یک روش جدید با اندازه گام $$h / 2$$ و جمله خطای $$O ( h ^ 4 )$$ است.

با فرایند مشابهی می‌توان دقت را بیش از این نیز افزایش داد. مجدداً $$h$$ را در $$I _ 1 (h/2)$$ با $$h / 2$$ جایگزین کرده و گام‌های مشابهی را طی می‌کنیم و خواهیم داشت:

$$\large \begin {equation} I [ f ] = \frac { 1 6 I _ 1 ( h / 4 ) - I _ 1 ( h / 2 ) }{ 1 5 } + O ( h ^ 6 ) = I _ 2 ( h / 4 ) + O ( h ^ 6 ) \end {equation} \;\;\; (44)$$

که در آن، یک روش جدید تعریف کردیم:

$$\large \begin {equation} I _ 2 ( h / 4 ) = \frac { 1 6 I _ 1 ( h / 4 ) - I _ 1 ( h / 2 ) }{ 1 5 } \end {equation} \;\;\; (45)$$

با اندازه گام $$h / 4$$ و جمله خطای $$O ( h ^ 6 )$$.

مثال انتگرال گیری با برون یابی ریچاردسون

انتگرال تابع گاوسی را با میانگین صفر و انحراف معیار واحد محاسبه می‌کنیم:

$$\large \begin {equation} I [ f ] = \int _ a ^ b f ( x ) \, d x = \frac { 1 }{ \sqrt { 2 \pi } } \int _ 0 ^ 3 e ^ { - x ^ 2 / 2 } \; d x \end {equation} \;\; (46)$$

در ابتدا، $$m = 0$$ و $$h_0=(b-a)/2^m=b-a=3$$ را قرار می‌دهیم و انتگرال را با قاعده ذوزنقه‌ای تقریب می‌زنیم:

$$\large \begin {equation} I [ f ] \approx I _ { 0 , 0 } = ( b - a ) \frac { f ( a ) + f ( b ) }{ 2 } = 0.6051 \end {equation} \;\; (47)$$

با $$n = 1$$ و $$h _ 1 = (b-a) / 2 = 1.5$$ و انتگرال را به صورت زیر تقریب می‌زنیم:‌

$$\large \begin {equation} I [ f ] \approx I _ { 1 , 0 } = \frac { h } { 2 }[ f ( a ) + 2 f ( a+ h ) + f ( b ) ] = 0 . 4 9 6 8 \end {equation} \;\; (48)$$

اکنون می‌توانیم یک روش جدید با برون یابی ریچاردسون برای تولید یک نتیجه دقیق‌تر به دست آوریم:

$$\large \begin {equation} I _ { 1 , 1 } = \frac { 4 I _ { 1 , 0 } - I _ { 0 , 0 } }{ 3 } = 0.4 6 0 7 2 \end {equation} \;\; (49)$$

این فرایند به صورت بازگشتی بیشتر برای تولید نتایج نشان داده شده در جدول زیر انجام می‌شود. هر وقت اندیس سطر $$m$$ به اندازه ۱ افزایش یابد، اندازه گام $$h_m=h/2^m$$ نصف می‌شود. و هر بار که اندیس ستون $$n$$ به یکی افزایش یابد، با استفاده از سطح بازگشتی بیشتر، روش جدید به دست می‌آید. تمام نتایج در ستون نهم با همان روش $$I_n$$ تولید می‌شوند، اما بر اساس اندازه‌های گام مختلف. همه نتایج سطر $$m$$اُم با روش‌های مختلف تولید می‌شوند، اما همه بر اساس همان اندازه گام $$h _ n = h / 2 ^ n$$ است. هر عنصر $$I _ {m,n}$$ در جدول بر اساس همسایه سمت چپ $$I_ {m ، n-1}$$ و همسایه سمت چپ بالای آن $$I_ {m-1 ، n-1}$$ توسط بازگشت زیر به دست می‌آید:

$$\large \begin {equation} I _ { m , n } = \frac { 4 ^ n I _ { m , n - 1 } - I _ { m - 1 , n - 1 } } { 4 ^ n - 1 } \end {equation} \;\; (50)$$

دقیق‌ترین نتیجه در هر سطح بازگشت روی قطر جدول قرار دارد. در این مورد خاص، $$I _ { 4 , 4 } = 0498650193$$ دقیق‌ترین نتیجه به دست آمده پس از $$n = 4$$ سطح بازگشت با استفاده از اندازه گام $$h_4=(b-a)/2^4=(b-a)/16$$ است.

سطر آخر خطاهای نسبی پنج روش $$I_{n,n},\;n=0,\cdots,4$$ در قطر را نشان می‌دهد.

کد متلب مورد استفاده برای تولید نتایج این جدول به صورت زیر است:

1    h=b-a;                     % initial step size
2    n=5;                       % number of recursions
3    m=n;
4    I=zeros(n);                % initializing array I
5    for i=1:m                  % for all rows
6        n=2^(i-1);
7        I(i,1)=(f(a)+f(b))/2;  
8        for k=1:n-1
9            I(i,1)=I(i,1)+f(a+k*h);
10        end
11        I(i,1)=h*I(i,1);       % composite trapezoidal rule
12        for j=2:i              % for all columns (Richardson extrapolation)
13            I(i,j)=(4^(j-1)*I(i,j-1)-I(i-1,j-1))/(4^(j-1)-1);
14        end
15        h=h/2;
16    end