درون یابی نیوتن — به زبان ساده (+ دانلود فیلم آموزش رایگان)
در آموزشهای قبلی مجله فرادرس، دیدیم که درونیابی لاگرانژ مبتنی بر $$ n + 1 $$ نقطه درونیابی $$ \{x_i,\,y_i=f(x_i),\;i=0,\cdots,n\} $$ است. این همه چیزی بود که برای محاسبه هر یک از چندجملهایهای اساسی یا پایه $$ l _ i ( x ) $$ نیاز داشتیم. اگر در صورت امکان بتوان از نقاط اضافهای استفاده کرد، باید همه چندجملهایهای پایه را در روش لاگرانژ تعویض کنیم. در مقایسه با روش لاگرانژ، در «درون یابی نیوتن» (Newton Interpolation)، وقتی نقاط بیشتری برای استفاده داشته باشیم، میتوانیم چندجملهایهای پایه بیشتر و ضرایب مربوط به آنها را محاسبه کنیم و چندجملهایهای پایه موجود و ضرایب آنها را بدون تغییر بگذاریم. بدین ترتیب، به دلیل وجود جملات اضافه، درجه چندجملهای درونیاب بزرگتر شده و ممکن است خطای تقریب کاهش پیدا کند (برای مثال، زمانی که چندجملهایهای مرتبه بالاتری را درونیابی کنیم).
درون یابی نیوتن
چندجملهایهای پایه درون یابی نیوتن به صورت زیر محاسبه میشوند:
$$ \large \begin {equation}
n _ 0 ( x ) = 1 , \; \; \; \; n _ i ( x ) = \prod _ { j = 0 } ^ { i - 1 } ( x - x _ j ) , \; \; \; \; \; \; ( i = 1 , \cdots , n )
\end {equation} $$
و چندجملهای درونیاب نیوتن به صورت زیر تشکیل میشود:
$$ \large \begin {eqnarray}
N _ n ( x ) & = & \sum _ { i = 0 } ^ n c _ i n _ i ( x ) = c _ 0 + \sum _{ i = 1 } ^ n c _ i \;
\left ( \prod _ { j = 0 } ^ { i - 1 } ( x - x _ j ) \right )
\nonumber \\
& = & c _ 0 + c _ 1 ( x - x _ 0 ) + c _ 2 ( x - x _ 0 ) ( x -x _ 1 ) + \cdots + c _ n \prod _ { j = 0 } ^ { n - 1 } ( x - x _ j )
\end {eqnarray} $$
لازم به ذکر است که آخرین نقطه $$ ( x _ n , y _ n ) $$ در هیچیک از چندجملهایهای پایه استفاده نمیشود، اما برای محاسبه آخرین ضریب $$ c _ n $$ به کار میرود. برای آنکه این چندجملهای $$ N _ n ( x ) $$ مرتبه $$n$$ از همه $$ n + 1 $$ نقطه $$ ( x _ i , y _ i ) $$ برای $$ i = 0 , \cdots , n $$ عبور کند، باید $$ n + 1 $$ معادله زیر برقرار باشند:
$$ \large \begin{equation}
\left \{ \begin {array} { l }
y _ 0 = N _ n ( x _ 0 ) = c _ 0 \\
y _ 1 = N _ n ( x _ 1 ) = c _ 0 + c _ 1 ( x _ 1 - x _ 0 ) \\
y _ 2 = N _ n ( x _ 2 ) = c_ 0 +c _ 1 (x _2 - x _ 0 ) +c _ 2 ( x _ 2 - x _ 0 ) (x _ 2 - x _ 1 ) \\
y _ 3 = N _ n ( x _ 3 ) = c _ 0 + c_ 1 ( x _ 3 - x _ 0 ) +c _ 2 (x _ 3 - x _ 0 ) ( x _ 3 - x _ 1 ) + c _ 3 ( x _ 3 - x _ 0 ) ( x _ 3 - x _ 1 ) ( x _ 3 - x _ 2 ) \\
\cdots \cdots \cdots \cdots \\
y _ n = N _ n ( x _ n ) = c _ 0 + c _ 1 ( x _ n - x _ 0 ) + c _ 2 ( x _ n - x _ 0 ) ( x _ n - x _ 1 ) + \cdots + c _ n \prod _ { i =0 } ^ { n - 1 } ( x _ n - x _ i )
\end {array} \right .
\end{equation} $$
که میتوان آن را به فرم ماتریسی زیر نوشت:
$$ \large \begin {equation}
\left [ \begin {array} { c c c c c c } 1 & & & & & \\ 1 & x _1 - x _ 0 & & & & \\
1 & x _ 2 - x _ 0 & ( x _ 2 - x _ 0 ) ( x _ 2 - x _ 1 ) & & & \\
1 & x _ 3 - x _ 0 & ( x _ 3 - x _ 0 ) ( x _ 3 - x _ 1 ) &( x _ 3 - x _ 0) ( x _ 3 - x _ 1 ) ( x _ 3 - x _ 2 ) & & \\
\vdots & \vdots & \vdots & \vdots & \ddots & \\
1 & x _ n - x _ 0 & ( x _ n - x_ 0 ) ( x _ n - x _ 1 ) & ( x _ n - x _ 0 ) ( x _ n - x _ 1 ) ( x _ n - x _ 2 ) & \cdots & \prod _ { i = 0 } ^ { n - 1 } ( x _n - x _ i )
\end {array} \right ]
\left [ \begin {array} { c } c _ 0 \\ c _ 1 \\ c _ 2 \\ c _ 3 \\ \vdots \\ c _ n \end {array} \right ]
= \left [ \begin {array} { c } y _ 0 \\ y _ 1 \\ y_ 2 \\ y _ 3 \\ \vdots \\ y _ n \end {array} \right ]
\end {equation} $$
تعداد $$ n + 1 $$ ضریب $$ c _ 0 $$، ... و $$ c _ n $$ را میتوان با حل این $$ n + 1 $$ معادله در دستگاه معادلات به دست آورد:
$$ \large \begin {eqnarray}
c _ 0 & = & y _ 0 = f ( x _ 0 ) = f [ x _ 0 ]
\nonumber \\
c _ 1 & = & \frac { y _ 1 - y _ 0 } { x _ 1 - x _ 0 } = f [ x _ 0 ,x _ 1 ]
\nonumber \\
c _ 2 & = & \frac { \frac { y _ 2 - y _ 0 } { x _ 2 - x _ 0 } - \frac { y _ 1 - y _ 0 } { x _ 1 - x _ 0 } } { x _ 2 - x _ 1 }
= f [ x _ 0 , x _ 1 , x _ 2 ]
\nonumber \\
& = & \frac { y _ 0 } { ( x _ 0 - x _ 2 ) ( x _ 0 - x _ 1 ) } + \frac { y _ 1 } { ( x _ 1 - x _ 0 ) ( x _ 1 - x _ 2 ) }
+ \frac { y _ 2 } { ( x _ 2 - x _ 0 ) ( x _ 2 - x _ 1 ) }
\nonumber \\
c _ 3 & = & \frac { \frac { \frac { y _ 3 - y _ 0 } { x _ 3 - x _ 0 } - \frac { y _ 1 - y _ 0 } { x _ 1 - x _ 0 } } { x _ 3 - x _ 1 } - \frac { \frac { y _ 2 - y _ 0 } { x _ 2 - x _ 0 } - \frac { y _ 1 - y _ 0 } { x _ 1 - x _ 0 } } { x _ 2 - x _ 1 } } { x _ 3 - x _ 2 }
= f [ x _ 0 , \cdots , x _ 3 ]
\nonumber \\
& = & \frac { y _ 0 } { ( x _ 0 - x _ 1 ) ( x _ 0 - x _ 2 ) ( x _ 0 - x _ 3 ) }
+ \frac { y _ 1 } { ( x _ 1 - x _ 0 ) ( x _1 - x _ 2 ) ( x _ 1 -x _ 3 ) }
\nonumber \\
& & + \frac { y _ 2 } { ( x _ 2 - x _ 0 ) ( x _ 2 - x _ 1 ) ( x_ 2 - x _ 3 ) }
+ \frac { y _ 3 } { ( x _ 3 - x _ 0 ) ( x _ 3 - x _ 1 ) ( x _ 3 -x _ 2 ) }
\nonumber \\
\cdots & & \cdots \cdots \cdots \cdots
\end {eqnarray} $$
در حالت کلی، داریم:
$$ \large \begin {equation}
c _ k = f [ x _ 0 , \cdots , x _ k ]
= \sum _ { j = 0 } ^ k \frac { f ( x _ j ) } { \prod _ { i = 0 , \; i \ne j } ^ k ( x _ j -x _ i ) } ,
\; \; \; \; \; ( k = 0 , \cdots , n )
\end {equation} $$
که فرم گسترش یافته از $$k$$اُمین تفاضل تقسیم شده $$ f [ x _ 0 , ... , x _ k ] $$ از $$k+1$$ نقطه اول است. اکنون درون یابی نیوتن را میتوان به صورت زیر نوشت:
$$ \large \begin {equation}
N _ n ( x ) = \sum _ { i = 0 } ^ n c _ i n _ i ( x ) = \sum _ { i = 0 } ^ n f [ x _ 0 , \cdots , x _ i ] n _ i ( x )
= f [ x _ 0 ] + \sum _ { i = 1 } ^ n \; f [ x _ 0 , \cdots , x _ i ] \left ( \prod _ { j = 0 } ^ { i - 1 } ( x - x _ j ) \right )
\end {equation} $$
به دلیل یکتایی درونیابی چندجملهای، این درون یابی نیوتن مانند درونیابیهای تابع توانی و لاگرانژ است: $$ V_ n ( x ) = L _ n ( x ) = P _ n ( x ) $$. اینها چندجملهایهای مرتبه $$n$$اُم مشابهی هستند، اما با چندجملهایهای پایه وزندار مختلفی با ضرایب متفاوت تشکیل میشوند.
چند نکته مهم در مورد درون یابی نیوتن به شرح زیر است:
نکته ۱: وقتی یک نقطه اضافه $$ ( x _ { n + 1 } , y _ { n + 1 } ) $$ داشته باشیم و از آن استفاده کنیم، همه چندجملهایهای پایه قبلی و ضرایب متناظر با آنها بدون تغییر باقی میمانند و تنها لازم است یک چندجملهای پایه جدید با درجه $$ n + 1 $$ به دست آوریم:
$$ \large \begin {equation}
n _ { n + 1 } ( x ) = n _ n ( x ) ( x - x _ n ) = \prod _ { i = 0 } ^ { n - 1 } ( x - x _ i ) ( x - x _ n )
= \prod _ { i = 0 } ^ n ( x -x _ i ) = \omega ( x )
\end {equation} $$
که دارای ضریب $$ (n+1)$$اُمین تفاضل تقسیم شده $$ c _ { n + 1 } = f [ x _ 0 , \cdots , x _ n , x _ { n + 1 } ] $$ است. چندجملهای درونیابی جدید با درجه $$ n + 1 $$ را میتوان با در نظر گفتن جمله اضافه در مجموع بالا به دست آورد:
$$ \large \begin {equation}
N _ { n + 1 } ( x ) = \sum _ { i = 0 } ^{ n + 1 } c _ i n _ i ( x )
= \sum _ { i = 0 } ^ n c _ i n _ i ( x ) + c _ { n + 1 } n_ { n + 1 } ( x )
= N _ n ( x ) + f [ x _ 0 , \cdots , x _ n , x _ { n + 1 } ] \omega ( x )
\end {equation} $$
از آنجا که $$ N _ { n + 1 } ( x ) $$ از نقطه جدید $$ ( x _ { n + 1 } , y _ { n + 1 } ) $$ میگذرد، داریم:
$$ \large \begin {equation}
f ( x _ { n + 1 } ) = N _ { n + 1 } ( x _ { n + 1 } )
= N _ n ( x _ { n + 1 } ) + f [ x _ 0 , \cdots , x _ n , x _ { n + 1 } ] \omega ( x )
\end {equation} $$
اما $$ x _ { n + 1 } $$ یک نقطه دلخواه است و میتوانیم آن را با $$ x $$ جایگزین کنیم. در این صورت، میتوان نوشت:
$$ \large \begin {equation}
f ( x ) = N _ n ( x ) + f [ x _ 0 , \cdots , x _ n , x ] \omega ( x ) \end {equation} $$
با مقایسه این تابع با تابع خطا، داریم:
$$ \large \begin {equation}
R _ n ( x ) = f ( x ) - N _ n ( x ) = \frac { f ^ { ( n + 1 ) } (\xi ) } { ( n + 1 ) ! } \omega ( x )
\end {equation} $$
مشاهده میکنیم که تابع خطا را به صورت زیر نیز میتوان نوشت:
$$ \large \begin {equation}
R _ n ( x ) = f [ x _ 0 , \cdots , x _ n , \, x ] \omega ( x )
\end {equation} $$
که در آن:
$$ \large \begin {equation}
f [ x _ 0 , \cdots , x _ n , x ] = \frac { f ^ { ( n + 1 ) } ( \xi ) } { ( n + 1 ) ! }
\end {equation} $$
نکته ۲: اگر همه $$ n + 1 $$ نقطه به یک موقعیت نزدیک باشند، یعنی $$ x _ i \to x _ 0 $$ ($$ i = 1, \cdots , n $$)، در حالت حدی، به نقطه مشابه $$ x _ 0 $$ تبدیل میشوند که $$ n $$ بار تکرار شده است. در نتیجه، داریم:
$$ \large \begin {equation}
f [ x _ 0 , \cdots , x _ n ] = \frac { f ^ { ( n ) } ( \xi ) } { n ! }
\; \; \; \stackrel { x _ i \rightarrow x _ 0 } { \Longrightarrow } \; \; \;
f [ x _ 0 , \cdots , x _ 0 ] = \frac { f ^ { ( n ) } ( x _ 0 ) } { n ! }
\end {equation} $$
و در نهایت، درون یابی نیوتن بر اساس $$ n + 1 $$ نقطه به شکل زیر در خواهد آمد:
$$ \large \begin {equation}
N _ n ( x ) = f [ x _ 0 ] + \sum _ { i = 1 } ^ n \left ( f [ x _ 0 , \cdots , x _ i ] \prod _ { j = 0 } ^ { i - 1 } ( x - x _ j ) \right )
\stackrel { x _ i \rightarrow x _ 0 } { \Longrightarrow }
f ( x _ 0 ) + \sum _ { i = 1 } ^ n \frac { f ^ { ( i ) } ( x _ 0 ) } { i ! } ( x -x _ 0 ) ^ i
\end {equation} $$
که $$ n + 1 $$ جمله اول بسط سری تیلور تابعی با خطای بریده شدن زیر است:
$$ \large \begin {equation}
R _ n ( x ) = \frac { f ^ { ( n + 1 ) } ( \xi ) } { ( n + 1 ) ! } \prod _ { i = 0 } ^ n ( x - x _ i )
\stackrel { x _ i \rightarrow x _ 0 } { \Longrightarrow }
\frac { f ^ { ( n + 1 ) } ( \xi ) } { (n + 1 ) ! } (x - x _ 0 ) ^{ n +1 }
\end {equation} $$
که در آن، $$ \xi$$ نقطهای بین $$ x $$ و $$ x _ 0 $$ است. میبینیم که سری تیلور در واقع حالت خاصی از درون یابی نیوتن در شرایط حدی است، وقتی که $$ n + 1 $$ نقطه به موقعیت مشابه $$ x _ 0 $$ نزدیک میشوند.
نکته ۳: اگر همه $$ n + 1 $$ نقطه $$ x _ 0 = a \le x _ 1 \le \cdots \le x _ { n - 1 } \le x _ n = b $$ فواصل برابری داشته باشند، یعنی داشته باشیم:
$$ \large \begin {equation}
x _ { i + 1 } - x _ i =h = \frac { b - a } { n } , \; \; \; \; \; ( i = 0 , \cdots , n - 1 )
\end {equation} $$
در نتیجه، درونیابی تفاضل تقسیم شده نیوتن را میتوان به فرم سادهتری نوشت. برای هر نقطه $$ x \in ( a , b ) $$، مقدار $$ c = ( x - x _ 0 ) / h $$ را به گونهای در نظر میگیریم که بتوان آن را به صورت $$ x = x _ 0 + c h $$ و $$ x - x _ i = ( x _ 0 + c h ) - ( x _ 0+i h ) = ( c - i ) h $$ نشان داد. اکنون چندجملهای نیوتن به صورت زیر نوشته میشود:
$$ \large \begin {eqnarray}
N ( x ) & = & N( x _ 0 + c h ) = f [ x _ 0 ] + \sum _ { i = 1 } ^ n f [ x _ 0 , \cdots , x _ i ] \prod _ { j = 0 } ^ { i - 1 } ( x - x _ j )
\nonumber \\
& = & f [ x _ 0 ] + f [ x _ 0 , x _ 1 ] c h + f [ x _ 0 , x _ 1 , x _ 2 ] c ( c - 1 ) h ^ 2
+ \cdots + f [ x _ 0 , \cdots , x _ n ] c ( c - 1 ) \cdots ( c -n + 1 ) h ^ n
\nonumber \\
& = & \sum _ { i = 0 } ^ n f [ x _ 0 , \cdots , x _ i ] c ( c - 1 ) ( c - 2 ) \cdots ( c - i + 1 ) h ^ i
= \sum _ { i = 0 } ^ n f [ x _ 0 , \cdots , x _ i ] \left ( \begin {array} { c } c \\ i \end {array} \right ) i ! \; h ^ i
\end {eqnarray} $$
که در آن:
$$ \large \begin {equation}
\left (
\begin {array} { c } c \\ i \end {array} \right )
= \frac { c ! } { i ! ( c - i ) ! } = \frac { c ( c - 1 ) ( c - 2 ) \cdots ( c - i + 1 ) } { i ! }
\end {equation} $$
مثال درون یابی نیوتن
تابع $$ y = f ( x ) = x \sin ( 2 x + \pi / 4 ) + 1 $$ را با یک چندجملهای درجه $$ n = 3 $$، بر اساس $$ n + 1 = 4 $$ نقطه زیر تقریب بزنید.
$$ \large \begin {equation}
\begin {array} { c | | c | c | c | c } \hline
i & 0 & 1 & 2 & 3 \\ \hline \hline
x _ i & - 1 & 0 & 1 & 2 \\ \hline
f ( x _i ) & 1.937 & 1.000 & 1.349 & -0.995 \\ \hline
\end {array}
\nonumber \\
\end {equation} $$
حل: بر اساس $$ f [ x _ i ] = f ( x _ i ), \; ( i = 0 , \cdots , n ) $$، میتوانیم همه تفاضلهای تقسیم شده دیگر را به صورت بازگشتی به فرم جدول زیر بنویسیم. در حالت کلی، $$ f [ x _ i , \cdots , x _ j ] $$ را میتوان بر اساس همسایه سمت چپ $$ f [ x _ { i + 1 } , \cdots , x _ j ] $$ و همسایه سمت چپ $$ f [ x _ i , \cdots , x _ { j - 1 } ] $$ به دست آورد:
$$ \large \begin {equation}
f [ x _ i , \cdots , x _ j ] = \frac { f [ x _ { i + 1 } , \cdots , x _ j ] - f [ x _ i , \cdots , x _ { j - 1 }] } { x _j - x _ i }
\nonumber \\
\end {equation} $$
$$ \large \begin {equation}
\begin {array} { c | | l |l | l | l } \hline
x _ i & 0 t h & 1 s t & 2 n d & 3 r d \\ \hline \hline
x _ 0 = - 1 & f [ x _ 0 ] = 1.937 & & & \\ \hline
x _ 1 = 0 & f [ x _ 1 ] = 1 . 0 0 0 & f [ x _ 0 , x _ 1 ] = - 0.937 & & \\ \hline
x _ 2 = 1 & f [ x _ 2 ] = 1.349 & f [ x _ 1 , x _ 2 ] = 0.349 & f [ x _ 0 ,x_1, x _ 2 ] = 0 .643 & \\ \hline
x _ 3 = 2 & f [ x _ 3 ] = - 0 . 9 9 5 & f [ x _ 2 , x _ 3 ] = - 2 . 3 43 & f [ x _ 1 , x_2, x _ 3 ] = - 1.3 4 6 & f [ x _ 0 , x_1 , x_2 , x _ 3 ] = - 0 . 6 6 3 \\ \hline
\end {array}
\nonumber \\
\end {equation} $$
ضرایب، چهار تفاضل تقسیم شده روی قطر هستند: $$ c_0=f[x_0]=f(x_0)=y_0=1.937$$، $$ c_1=f[x_0,x_1]=-0.937$$، $$ c_2=f[x_0,x_1,x_2]=0.643$$ و $$ c_3=f[x_0,x_1,x_2,x_3]=-0.663 $$. یک راه دیگر برای نشان دادن ضرایب، استفاده از فرم گسترش یافته است:
$$ \large \begin {equation}
c _ i = f [ x _ 0 , \cdots , x _ i ] = \sum _ { j = 0 } ^ n \frac { f ( x _ j ) } { \prod _ { i = 0 , \, i \ne j } ^ n ( x _ j - x _ i ) } \end {equation} $$
در نهایت، چندجملهای درونیاب نیوتن به صورت زیر به دست میآید:
$$ \large \begin {eqnarray}
N _ 3 ( x ) & = & \sum _ { i = 0 } ^ 3 c _ i n _ i ( x )
\nonumber \\
& = & 1.937 - 0.937 ( x + 1 ) + 0.643 ( x + 1 ) ( x - 0 ) - 0.663 ( x + 1 ) ( x - 0 ) ( x - 1 )
\nonumber \\
& = & 1 + 0.369 \, x + 0.643 \, x ^ 2 - 0 . 663 \, x ^ 3
\nonumber
\end {eqnarray} $$
درونیابیهای چندجملهای با استفاده از روش سری توانی، لاگرانژ و نیوتن دقیقاً با هم برابرند: $$ P _ 3 ( x ) = L _ 3 ( x ) = N _ 3 ( x )$$. این مورد یکتایی درونیابی چندجملهای را نشان میدهد که در شکلهای زیر برای $$ y = f ( x ) $$ نشان داده شده است.
پیادهسازی درون یابی نیوتن در متلب
کد متلب روش دون یابی نیوتن در ادامه ارائه شده است. ضرایب $$ c _ i = f [ x _ 0 , \cdots , x _ n ] , \; ( i = 0 , \cdots , n ) $$ را میتوان با استفاده از فرم گسترش یافته یا فرم جدول به صورت بازگشتی تولید کرد.
1function [v N]=NI(u,x,y) % Newton's Interpolation
2 % vectors x and y contain n+1 points and the corresponding function values
3 % vector u contains all discrete samples of the continuous argument of f(x)
4 n=length(x); % number of interpolating points
5 k=length(u); % number of discrete sample points
6 v=zeros(1,k); % Newton interpolation
7 N=ones(n,k); % all n Newton's polynomials (each of m elements)
8 N(1,:)=y(1); % first Newton's polynomial
9 v=v+N(1,:);
10 for i=2:n % generate remaining Newton's polynomials
11 for j=1:i-1
12 N(i,:)=N(i,:).*(u-x(j));
13 end
14 c=DividedDifference(x,y,i) % get the ith coefficient c_i
15 v=v+c*N(i,:); % weighted sum of all Newton's polynomials
16 end
17end
18
19function dd=DividedDifference(x,y,i) % generate f[x_0,...,x_i] in expanded form
20 dd=0;
21 for k=1:i % loop for summation
22 de=1;
23 for l=1:i % loop for product
24 if k~=l
25 de=de*(x(k)-x(l));
26 end
27 end
28 dd=dd+y(k)/de; % ith coefficient c_i
29 end
30end
31
32function dd=DividedDifferenceMatrix(x,y) % generate divided difference matrix
33 n=length(x); % the coefficients are along diagonal
34 dd=zeros(n); % matrix of divided differences
35 dd(:,1)=y;
36 for i=1:n
37 fprintf('%6.3f\t',dd(i,1))
38 for j=2:i
39 dd(i,j)=(dd(i,j-1)-dd(i-1,j-1))/(x(i)-x(i-j+1));
40 fprintf('%6.3f\t',dd(i,j));
41 end
42 fprintf('\n');
43 end
44end
اگر این مطلب برایتان مفید بوده است، آموزشهای زیر نیز به شما پیشنهاد میشوند:
- مجموعه آموزشهای محاسبات عددی
- آموزش محاسبات عددی (مرور و حل مساله)
- مجموعه آموزشهای دروس ریاضیات
- آموزش محاسبات عددی با MATLAB
- درون یابی اسپلاین — از صفر تا صد
- درون یابی هرمیت — از صفر تا صد
- درون یابی خطی — به زبان ساده
^^