Matlab Tran Van Chinh Chuong 4 Noi Suy Va Xap Xi Ham - Tài Liệu Text

Có thể bạn quan tâm

Tải bản đầy đủ (.pdf) (31 trang)

Trang chủ

Giáo án - Bài giảng

Cao đẳng - Đại học

matlab tran van chinh chuong 4 noi suy va xap xi ham

Bạn đang xem bản rút gọn của tài liệu. Xem và tải ngay bản đầy đủ của tài liệu tại đây (410.04 KB, 31 trang )

CHƯƠNG 3: NỘI SUY VÀ XẤP XỈ HÀM§1. NỘI SUY LAGRANGE Trong thực tế nhiều khi ta cần tính giá trị của hàm y = f(x) tại một giá trị x trong một đoạn [a, b] nào đó mà chỉ biết một số nhất định các giá trị của hàm tại một số điểm cho trước. Các giá trị này được cung cấp qua thực nghiệm hay tính tốn. Vì vậy nảy sinh vấn đề tốn học là trên đoạn a ≤ x ≤ b cho một loạt các điểm xi ( i = 0, 1, 2...) và tại các điểm xi này giá trị của hàm là yi = f(xi) đã biết và ta cần tìm y = f(x) dựa trên các giá trị đã biết đó. Lúc đó ta cần tìm đa thức : Pn(x) = aoxn + a1xn‐1 + …+an‐1x + an sao cho Pn(xi) = f(xi) = yi. Đa thức Pn(x) được gọi là đa thức nội suy của hàm y=f(x). Ta chọn đa thức để nội suy hàm y = f(x) vì đa thức là loại hàm đơn giản, ln có đạo hàm và ngun hàm. Việc tính giá trị của nó theo thuật tốn Horner cũng đơn giản. Bây giờ ta xây dựng đa thức nội suy kiểu Lagrange. Gọi Li là đa thức: ( x − x0 )...( x − xi −1 )( x − xi + 1 )...( x − x n )Li = ( xi − x 0 )...( xi − xi −1 )( x i − x i + 1 )...( x i − x n ) Rõ ràng là Li(x) là một đa thức bậc n và : j=i⎧1L i (x j ) = ⎨ 0j≠i⎩Ta gọi đa thức này là đa thức Lagrange cơ bản. Bây giờ ta xét biểu thức : n Pn ( x) = ∑ f( x i )L i ( x) i =0Ta thấy Pn(x) là một đa thức bậc n vì các Li(x) là các đa thức bậc n và thoả mãn điều kiện Pn(xi) = f(xi) = yi. Ta gọi nó là đa thức nội suy Lagrange. Với n = 1 ta có bảng x x0 x1 y y0 y 1 Đa thức nội suy sẽ là : P1(x) = yoL0(x) + y1L1(x1) x − x0x − x1L0 = L 1 = x 0 − x1x1 − x 0210CuuDuongThanCong.com /> nên P1 ( x) = y 0x − x0x − x1+ y1 x1 − x 0x 0 − x1Như vậy P1(x) là một đa thức bậc nhất đối với x Với n = 2 ta có bảng x x0 x1 x2 y y 0 y1 y2 Đa thức nội suy sẽ là : P2(x) = yoL0(x) + y1L1(x1) + y2L2(x2) ( x − x1 )( x − x 2 ) L0 = ( x 0 − x1 )( x0 − x 2 )( x − x0 )( x − x 2 )L1 = ( x1 − x0 )( x1 − x 2 )( x − x 0 )( x − x1 ) L2 =( x 2 − x 0 )( x 2 − x1 )Như vậy P1(x) là một đa thức bậc hai đối với x. Ta xây dựng hàm lagrange() để thực hiện việc nội suy hàm theo thuật toán Lagrange: function [l, L] = lagrange(x, y) %Dua vao : x = [x0 x1 ... xn], y = [y0 y1 ... yn] %ket qua: l = He so cua da thuc Lagrange bac n % L = Da thuc Lagrange n = length(x) ‐ 1; %bac cua da thucl l = 0; for m = 1:n + 1 p = 1; for k = 1:n + 1 if k ~= m p = conv(p, [1 ‐x(k)])/(x(m) ‐ x(k)); end end L(m, :) = p; %da thuc Lagrange l = l + y(m)*p; end 211CuuDuongThanCong.com /> Cho hàm dưới dạng bảng: x ‐2 ‐1 1 2 y ‐6 0 0 6 và tìm y(2.5) ta dùng chương trình ctlagrange.m: clear all, clc x = [‐2 ‐1 1 2]; y = [‐6 0 0 6]; l = lagrange(x, y); yx = polyval(l, 2.5) §2. NỘI SUY NEWTON Bây giờ ta xét một cách khác để xây dựng đa thức nội suy gọi là phương pháp Newton. Trước hết ta đưa vào một khái niệm mới là tỉ hiệu Giả sử hàm y = y(x) có giá trị cho trong bảng sau: x x0 x1 x2 … xn‐1 xn y y0 y1 y2 … yn‐1 yn Tỉ hiệu cấp 1 của y tại xi, xj là : yi − y j y[x i , x j ] = xi − x j Tỉ hiệu cấp hai của y tại xi, xj, xk là : y[x i , x j ] − y[x j , x k ]y[xi , x j , x k ] = xi − xkv.v. Với y(x) = Pn(x) là một đa thức bậc n thì tỉ hiệu cấp 1 tại x, x0 : P ( x) − Pn ( x0 ) Pn [x , x0 ] = nx − x0là một đa thức bậc (n ‐ 1). Tỉ hiệu cấp 2 tại x, x0, x1 : P [x , x0 ] − Pn [x0 , x1 ]Pn [x , x 0 , x1 ] = n x − x1là một đa thức bậc (n‐2) v.v và tới tỉ hiệu cấp (n + 1) thì : 212CuuDuongThanCong.com /> Pn[ x, xo,.., xn] = 0 Từ các định nghĩa tỉ hiệu ta suy ra : Pn(x) = Pn(x0) + ( x‐ x0)Pn[x, xo] Pn[x, x0] = Pn[x0, x1] + ( x ‐ x1)Pn[x, xo,x1] Pn[x, xo, x1] = Pn[x0, x1, x2] + ( x ‐ x2)Pn[x, xo, x1, x2] ............ Pn[x, xo,.., xn‐1] = Pn[x0, x1,.., xn] + ( x ‐ xn)Pn[x, xo,.., xn] Do Pn[ x, xo,.., xn] = 0 nên từ đó ta có : Pn(x) = Pn(x0) + (x ‐ x0)Pn[xo, x1] + (x ‐ x0)(x ‐ x1)Pn[x0, x1, x2] +… +(x ‐ x0)…(x ‐ xn‐1)Pn[x0,…, xn] Nếu Pn(x) là đa thức nội suy của hàm y = f(x) thì: Pn(xi) = f(xi) = yi với i = 0 ÷ n Do đó các tỉ hiệu từ cấp 1 đến cấp n của Pn và của y là trùng nhau và như vậy ta có : Pn(x) = y0 + (x ‐ x0)y[x0, x1] + (x ‐ x0)(x ‐ x1)y[x0, x1, x2] + .. + (x ‐ x0)(x ‐ x1)...(x ‐ xn‐1)y[x0,..,xn] Đa thức này gọi là đa thức nội suy Newton tiến xuất phát từ nút x0 của hàm y = f(x). Ngồi đa thức tiến cịn có đa thức nội suy Newton lùi xuất phát từ điểm xn có dạng như sau : Pn(x) = yn + (x ‐ xn)y[xn, xn‐1] + (x ‐ xn)(x ‐ xn‐1)y[xn, xn‐1,xn‐2] +..+ (x ‐ xn)(x ‐ xn‐1)...(x ‐ x1)y[xn,.., x0] Trường hợp các nút cách đều thì xi = x0 + ih với i = 0, 1,.., n. Ta gọi sai phân tiến cấp 1 tại i là : ∆yi = yi+1 ‐ yi và sai phân tiến cấp hai tại i: ∆2yi = ∆(∆yi) = yi+2 ‐ 2yi+1 + yi ......... và sai phân tiến cấp n là : ∆nyi = ∆(∆n‐1yi) Khi đó ta có: ∆y 0 y[x 0 , x 1 ] = h∆2 y 0 y[x 0 , x 1 , x 2 ] =2h 2 ........... 213CuuDuongThanCong.com /> ∆n y 0 y[x 0 , x 1 , x 2 ,..., x n ] = n! h nBây giờ đặt x = x0 + ht trong đa thức Newton tiến ta được: t( t − 1) 2t( t − 1) ⋅ ⋅ ⋅ ( t − n + 1) n Pn ( x 0 + ht) = y 0 + t∆y 0 +∆ y0 + ⋅ ⋅ ⋅ +∆ y0 2!n!thì ta nhận được đa thức Newton tiến xuất phát từ x0 trong trường hợp nút cách đều. Với n = 1 ta có : P1(x0 + ht) = y0 + ∆y0 Với n = 2 ta có: t( t − 1) 2Pn ( x 0 + ht) = y 0 + t∆y 0 +∆ y0 2!Một cách tương tự ta có khái niệm các sai phân lùi tại i: ∇yi = yi ‐ yi‐1 ∇2yi = ∇(∇yi) = yi ‐ 2yi‐1 + yi‐2 ......... ∇nyi = ∇(∇n‐1yi) và đa thức nội suy Newton lùi khi các điểm nội suy cách đều: t( t + 1) 2t( t + 1) ⋅ ⋅ ⋅ ( t + n − 1) nPn ( x 0 + ht) = y n + t∇y n +∇ yn + ⋅ ⋅ ⋅ +∇ yn 2!n!Ta xây dựng hàm newton() để nội suy: function [n,DD] = newton(x,y) %Dua vao : x = [x0 x1 ... xN] % y = [y0 y1 ... yN] %Lay ra: n = he so cua da thuc Newton bac N N = length(x) ‐ 1; DD = zeros(N + 1, N + 1); DD(1:N + 1, 1) = yʹ; for k = 2:N + 1 for m = 1: N + 2 ‐ k DD(m,k) = (DD(m + 1, k ‐ 1) ‐ DD(m, k ‐ 1))/(x(m + k ‐ 1) ‐ x(m)); end end a = DD(1, :); n = a(N+1); for k = N:‐1:1 214CuuDuongThanCong.com /> n = [n a(k)] ‐ [0 n*x(k)]; end Cho hàm dưới dạng bảng: x ‐2 ‐1 1 2 4 y ‐6 0 0 6 60 Ta dùng chương trình ctnewton.m để nội suy: clear all, clc x = [‐2 ‐1 1 2 4]; y = [‐6 0 0 6 60]; a = newton(x, y) yx = polyval(a, 2.5) §3. NỘI SUY AITKEN ‐ NEVILLE Một dạng khác của đa thức nội suy được xác định bằng thuật tốn Aitken ‐ Neville. Giả sử ta có n điểm đã cho của hàm f(x). Như vậy qua hai điểm x0 và x1 ta có đa thức nội suy Lagrange của hàm f(x) được viết dưới dạng: y0 x0 − xy x1 − xP01 ( x) = 1 x1 − x 0Đây là một đa thức bậc 1: x − x1x − x0P01 ( x) = y 0 + y1x 0 − x1x1 − x 0Khi x = x0 thì: P01 ( x 0 ) =y0y1x0 − x0x1 − x 0= y0 x1 − x 0Khi x = x1 thì: y 0 x 0 − x1y x1 − x1P01 ( x1 ) = 1= y1 x1 − x 0Đa thức nội suy Lagrange của f(x) qua 3 điểm x0, x1, x2 có dạng: 215CuuDuongThanCong.com /> P01 ( x) x0 − xP ( x) x 2 − xP012 ( x) = 12 x2 − x0và là một đa thức bậc 2: ( x − x1 )( x − x 2 )( x − x 0 )( x − x 2 )( x − x 0 )( x − x1 ) + y2+ y1 P012 ( x) = y 0( x 0 − x1 )( x 0 − x 2 )( x1 − x 0 )( x1 − x 2 )( x 2 − x 0 )( x 2 − x1 )Khi x = x0 thì: y0x0 − x0P ( x) x 2 − x 0P012 ( x0 ) = 12= y0 x 2 − x0Khi x = x1 thì: y 1 x 0 − x1y x 2 − x1P012 ( x1 ) = 1= y1 x2 − x0Khi x = x2 thì: P01 ( x 2 ) x0 − x 2y2x2 − x2P012 ( x 2 ) == y2 x2 − x0Tổng quát đa thức nội suy Lagrange qua n điểm là: P01..( n −1) ( x) x 0 − xP12..n ( x) x n − x P012..n ( x) =x2 − x0Như vậy ta có thể dùng phép lặp để xác định lần lượt các đa thức Lagrange. Sơ đồ tính tốn như vậy gọi là sơ đồ Neville ‐ Aitken. Ta xây dựng hàm aitkenneville() để nội suy: function a = aitkenneville(xData, yData, x) % Tra ve gia tri noi suy tai x. % Cu phap: y = aitkenneville(xData, yData, x) n = length(xData); y = yData; for k = 1:n‐1 y(1:n‐k) = ((x ‐ xData(k+1:n)).*y(1:n‐k)... + (xData(1:n‐k) ‐ x).*y(2:n‐k+1))... ./(xData(1:n‐k) ‐ xData(k+1:n)); 216CuuDuongThanCong.com /> end a = y(1); Cho các cặp số (1, 3), (2, 5), (3, 7), (4, 9) và (5, 11), để tìm y tại x = 2.5 ta dùng chương trình ctaitkennevile.m: clear all, clc x = [1 2 3 4]; y = [3 5 7 9]; yx = aitkenneville(x, y, 2.5) §4. NỘI SUY BẰNG ĐƯỜNG CONG SPLINE BẬC BA Khi số điểm cho trước dùng khi nội suy tăng, đa thức nội suy có dạng sóng và sai số tăng. Ta xét hàm thực: 1 f31(x) = 1 + 8x 2và nội suy nó bằng thuật tốn Newton nhờ chương trình cttestintp.m %Noi suy Newton x1 = [‐1 ‐0.5 0 0.5 1.0]; y1 = f31(x1); n1 = newton(x1,y1) x2 = [‐1 ‐0.75 ‐0.5 ‐0.25 0 0.25 0.5 0.75 1.0]; y2 = f31(x2); n2 = newton(x2,y2) x3 = [‐1 ‐0.8 ‐0.6 ‐0.4 ‐0.2 0 0.2 0.4 0.6 0.8 1.0]; y3 = f31(x3); n3 = newton(x3,y3) xx = [‐1:0.02: 1]; %pham vi noi suy yy = f31(xx); %ham thuc yy1 = polyval(n1, xx); %ham xap xi qua 5 diem yy2 = polyval(n2, xx); %ham xap xi qua 9 diem yy3 = polyval(n3, xx); %ham xap xi qua 11 diem subplot(221) plot(xx, yy, ʹk‐ʹ, xx, yy1, ʹbʹ) subplot(224) 217CuuDuongThanCong.com /> plot(xx, yy1‐yy, ʹrʹ, xx, yy2‐yy, ʹgʹ, xx, yy3‐yy,ʹbʹ) %do thi sai so subplot(222) plot(xx,yy,ʹk‐ʹ, xx, yy2, ʹbʹ) subplot(223) plot(xx, yy, ʹk‐ʹ, xx, yy3, ʹbʹ) yvà nhận được kết quả. fi,i+1 fi‐1,i Để tránh hiện tượng sai số lớn khi số điểm mốc tăng ta dùng nội suy nối trơn(spline). Trên các đoạn nội suy ta yi‐1 yi yi+1 thay hàm bằng một đường cong. Các xđường cong này được ghép trơn tại các xi xi+1 xi‐1 điểm nối. Ta chọn các đường cong này là hàm bậc 3 vì hàm bậc 1 và bậc hai khó bảo đảm điều kiện nối trơn. Cho một loạt giá trị nội suy (x1, y1),…,(xi, yi),…,(xn, yn). Trên mỗi đoạn ta có một hàm bậc 3. Như vậy giữa nút i và (i +1) ta có hàm fi,i+1(x), nghĩa là ta dùng (n ‐ 1) hàm bậc 3 f1,2(x), f2,3(x),…, fn‐1,n(x) để thay thế cho hàm thực. Hàm fi,i+1(x) có dạng: (1) fi,i+1(x) = ai + bi(x ‐ xi) + ci(x ‐ xi)2 + di(x ‐ xi)3 Hàm này thoả mãn: (3) fi,i+1(xi) = ai = yi 32 fi ,i+1 (xi+1 ) = di h i + c i h i + bi h i + a i = y i+1 (4) fi′,i+1 (xi ) = bi (5) fi′,i+1 (xi+1 ) = 3di h i2 + 2c i h i + bi fi′′,i+1 (x i ) = 2c i = y′′i fi′′,i+1 (xi+1 ) = 6di h i + 2c i = y′′i+1 (6) (7) (8) Muốn nối trơn ta cần có đạo hàm bậc nhất liên tục và do đó: fi′′−1,i (x i ) = fi′′,i+1 (x i ) = k i Lúc này các giá trị k chưa biết, ngoại trừ k1 = kn = 0(ta các các mút là điểm uốn). Điểm xuất phát để tính các hệ số của fi,i+1(x) là biểu thức của fi′′,i+1 (xi ) . Sử dụng nội suy Lagrange cho hai điểm ta có: fi′′,i+1 (x i ) = k i L i (x) + k i+1L i+1 (x) Trong đó: 218CuuDuongThanCong.com /> Li (x) =x − x i +1x i − x i +1Li+1 (x) =x − xix i +1 − x i Do vậy: k i (x − x i+1 ) − k i+1 (x − x i ) x i − x i +1Tích phân biểu thức trên hai lần theo x ta có: k i (x − xi+1 )3 − k i+1 (x − xi )3 fi ,i+1 (xi ) =+ A(x − xi+1 ) − B(x − xi ) 6(xi − xi+1 )Trong đó A và B là các hằng số tích phân Số hạng cuối trong phương trình trên thường được viết là Cx + D. Đặt C = A ‐ B và D = ‐Axi+1 + Bxi để dễ dàng tính tốn. Từ điều kiện fi,i+1(xi) = yi ta có: k i (xi − xi+1 )3+ A(xi − xi+1 ) = y i 6(xi − x i+1 )nên: yik (x − x i+1 )A=− i i x i − x i +16 fi′′,i+1 (x i ) =Tương tự, điều kiện fi,i+1(xi+1) = yi+1 cho ta: y i +1k (x − xi+1 )B=− i +1 i x i − x i +16Kết quả là: ⎤k ⎡ (x − xi+1 )3− (x − xi+1 )(xi − xi+1 ) ⎥fi ,i+1 (xi ) = i ⎢6 ⎣ x i − x i +1⎦⎤k i+1 ⎡ (x − xi )3−−−−(xx)(xx)iii +1 ⎥ 6 ⎢⎣ xi − xi+1⎦ y i (x − xi+1 ) − y i+1 (x − xi )x i − x i +1Đạo hàm cấp 2 ki tại các nút bên trong được tính từ điều kiện: fi′−1,i (x i ) = fi′,i+1 (x i ) +Sau khi biến đổi ta có phương trình: k i−1 (xi−1 − xi ) + 2k i (xi−1 − xi+1 ) + k i+1 (xi − xi+1 ) ⎛ y − y i y i − y i +1 ⎞= 6 ⎜ i −1−⎟⎝ x i −1 − x i x i − x i + 1 ⎠Khi các điểm chia cách đều (xi+1 ‐ xi) = h ta có: 219CuuDuongThanCong.com /> 6i = 2, 3,…, n ‐ 1 ( yi−1 − 2yi + yi+1 ) h2Ta xây dựng hàm cubicspline() để nội suy: function y = cubicspline(xData, yData, x) %Ham nay xap xi bang da thuc bac 3 spline %Cu phap: [yi,f] = cubicspline(xData, yData, x) n = length(xData); c = zeros(n‐1, 1); d = ones(n, 1); e = zeros(n‐1, 1); k = zeros(n, 1); c(1:n‐2) = xData(1:n‐2) ‐ xData(2:n‐1); d(2:n‐1) = 2*(xData(1:n‐2) ‐ xData(3:n)); e(2:n‐1) = xData(2:n‐1) ‐ xData(3:n); k(2:n‐1) = 6*(yData(1:n‐2) ‐ yData(2:n‐1))... ./(xData(1:n‐2) ‐ xData(2:n‐1))... ‐ 6*(yData(2:n‐1) ‐ yData(3:n))... ./(xData(2:n‐1) ‐ xData(3:n)); [c, d, e] = band3(c, d e); k = band3sol(c, d, e, k); i = findseg(xData, x); h = xData(i) ‐ xData(i+1); y = ((x ‐ xData(i+1))^3/h ‐ (x ‐ xData(i+1))*h)*k(i)/6.0... ‐ ((x ‐ xData(i))^3/h ‐ (x ‐ xData(i))*h)*k(i+1)/6.0... + yData(i)*(x ‐ xData(i+1))/h... ‐ yData(i+1)*(x ‐ xData(i))/h; Ta có chương trình ctcubicspline.m dùng nội suy: clear all, clc x1 = 0:0.1:5; y1 = (x1+1).^2; while 1 x = input(ʹx = ʹ); if isempty(x) fprintf(ʹKet thucʹ); break k i−1 + 4k i + k i+1 =220CuuDuongThanCong.com /> end y = cubicspline(xData, yData, x) fprintf(ʹ\nʹ) end §5. NỘI SUY BẰNG ĐA THỨC CHEBYSHEV Khi nội suy bằng đa thức Newton hay Lagrange, nghĩa là thay hàm thực bằng đa thức xấp xỉ, có khoảng cách cách đều thì sai số giữa đa thức nội suy và hàm thực có xu hướng tăng tại hai mút nội suy. Ta thấy rõ điều này khi chạy chương trình cttestintp.m. Do vậy ta nên chọn các điểm mốc nội suy ở hai mút dày hơn ở giữa. Một trong những cách chọn phân bố các điểm mốc là hình chiếu lên trục x của các điểm cách đều trên đường tròn 1tâm tại điểm giữa của đoạn nội suy. Như vậy với ‐1 x′1đoạn nội suy [‐1, 1] ta có: 2n + 1 − 2k x′k = cosπ k = 1, 2,…,n (1) 2(n + 1)Với đoạn nội suy [a, b] bất kì: b−ab+a b−a2n + 1 − 2ka+bxk =x′k +cos =π+ k = 1, 2,…,n (2) 2222(n + 1)2Các nút nội suy này được gọi là các nút Chebyshev. Đa thức nội suy dựa trên các nút Chebyschev gọi là đa thức nội suy Chebyshev. Ta xét hàm thực: 1f(x) = 1 + 8x 2Ta chọn số nút nội suy lần lượt là 5, 9, 11 và xây dựng các đa thức Newton (hay Lagrange) c4(x), c8(x) và c10(x) đi qua các nút này và vẽ đồ thị của hàm thực cũng như sai số khi nội suy bằng chương trình ctcomchebynew.m với các N khác nhau. x1 = [‐1 ‐0.5 0 0.5 1.0]; y1 = f31(x1); n1 = newton(x1,y1); xx = [‐1:0.02: 1]; %pham vi noi suy yy1 = polyval(n1,xx); %ham xap xi qua 5 diem yy = f31(xx); %ham thuc 221CuuDuongThanCong.com /> subplot(221) plot(xx,yy,ʹk‐ʹ, x, y, ʹoʹ, xx, yy1, ʹbʹ); title(ʹNewtonʹ) subplot(223) plot(xx, yy1‐yy, ʹrʹ) %do thi sai so N = 4; k = [0:N]; x = cos((2*N + 1 ‐ 2*k)*pi/2/(N + 1)); y = f31(x); c = newton(x, y) %da thuc noi suy dua tren cac nut Chebyshev xx = [‐1:0.02: 1]; %doan noi suy yy = f31(xx); %do thi ham thuc yy1 = polyval(c, xx); %do thi ham xap xi subplot(222) plot(xx, yy, ʹk‐ʹ, x, y, ʹoʹ, xx, yy1, ʹbʹ) title(ʹChebyshevʹ) subplot(224) plot(xx, yy1‐yy, ʹrʹ) %do thi sai so Khi tăng số điểm mốc, nghĩa là tăng bậc của đa thức Chebyschev, sai số giảm. Đa thức Chebyshev bậc n được xác định bằng: Tn+1(xʹ) = cos((n+1)arccos(xʹ)) (3) và các nút Chebyshev cho bởi (1) là nghiệm của (3). Ta có: Tn +1 (x′) = cos(arccos(x′) + narccos(x′))= cos(arccos(x′))cos(narccos(x′) − sin(arccos(x′))sin(narccos(x′)) = x′T n(x′) + 0.5 ⎣⎡cos((n + 1)arccos(x′) − cos((n − 1)arccos(x′)⎦⎤= x′T n(x′) + 0.5T n +1(x′) − 0.5T n −1(x′)nên: Tn +1(x′) = 2xT n(x′) − T n −1(x′) n ≥ 1 và T0(xʹ) = 1 T1(xʹ) = cos(arccos(xʹ) = xʹ Các đa thức Chebyshev đến bậc 6 là: T0(x) = 1 T1(xʹ) = xʹ T2(xʹ) = 2xʹ2 ‐ 1 (4) (5) 222CuuDuongThanCong.com /> T3(xʹ) = 4xʹ3 ‐ 3xʹ T4(xʹ) = 8xʹ4 ‐ 8ʹx2 + 1 T5(xʹ) = 16xʹ5 ‐ 20ʹx3 + 5xʹ T6(xʹ) = 32xʹ6 ‐ 48xʹ4 + 18xʹ2 ‐ 1 T7(xʹ) = 64xʹ7 ‐ 112xʹ5 + 56xʹ3 ‐ 7xʹ Hàm f(x) được xấp xỉ bằng: N f(x) = ∑ d m Tm (x′) x′=m =02 ⎛ a+b ⎞⎜ x−⎟b −a ⎝2 ⎠ (6) Trong đó: 1 n1 n d0 =f(x k ) ∑ f(xk )T0 (x′k ) = n + 1 ∑n + 1 k =0k =02 ndm =∑ f(xk )Tm (x′k )n + 1 k =0 2 nm(2n + 1 − 2k)=∑ f(xk )cos 2(n + 1) π m = 1,2,...,nn + 1 k =0Ta xây dựng hàm cheby() để tìm đa thức nội suy Chebyshev: function [c, x, y] = cheby(f, N, a, b) %vao : f = ten ham tren doan [a, b] %Ra: c = Cac he so cua da thuc Newton bac N % (x,y) = cac nut Chebyshev if nargin == 2 a = ‐1; b = 1; end k = [0: N]; theta = (2*N + 1 ‐ 2*k)*pi/(2*N + 2); xn = cos(theta); %pt.(1) x = (b ‐ a)/2*xn +(a + b)/2; %pt.(2) y = feval(f,x); d(1) = y*ones(N + 1,1)/(N+1); for m = 2: N + 1 cos_mth = cos((m‐1)*theta); d(m) = y*cos_mthʹ*2/(N + 1); %pt.(7) end xn = [2 ‐(a + b)]/(b ‐ a); %nghich dao cua t. (2) (7) (8) 223CuuDuongThanCong.com /> T_0 = 1; T_1 = xn; %pt.(5) c = d(1)*[0 T_0] +d(2)*T_1; %pt.(6) for m = 3: N + 1 tmp = T_1; T_1 = 2*conv(xn,T_1) ‐[0 0 T_0]; %pt.(4) T_0 = tmp; c = [0 c] + d(m)*T_1; %pt.(6) end Để tìm đa thức Chebyshev dùng xấp xỉ hàm f(x) =1 ta dùng chương 1 + 8x 2trình ctcheby.m: clear all, clc N = 2; a = ‐2; b = 2; [c, x1, y1] = cheby(ʹf31ʹ, N, a, b) %da thuc Chebyshev %so sanh voi da thuc Lagrange/Newton k = [0:N]; xn = cos((2*N + 1 ‐ 2*k)*pi/2/(N + 1));%pt.(1):nut Chebyshev x = ((b‐a)*xn +a + b)/2; %pt.(2) y = f31(x); n = newton(x, y) l = lagrange(x, y) §6. XẤP XỈ HÀM BẰNG PHÂN THỨC HỮU TỈ Xấp xỉ Padé dùng để xấp xỉ hàm f(x) tại x0 bằng hàm hữu tỉ: Q (x − x 0 )Pm ,n (x − x 0 ) = m D n (x − x 0 )q 0 + q 1 (x − x0 ) + q 2 (x − x0 )2 + L + q m (x − x0 )m = 1 + d1 (x − x0 ) + d 2 (x − x0 )2 + L + d n (x − x0 )nvới m = n hay m = n + 1 Trong đó f(x0), fʹ(x0),…, f(m+n)(x0) đã cho Trước hết ta khai triển Taylor hàm f(x) tại x = x0 đến bậc (m + n). (1) 224CuuDuongThanCong.com /> f(x) ≈ Tm + n (x) = f(x0 ) + f ′(x0 )(x − x0 )f′′(x0 )f (m + n) (x0 )2+(x − x0 ) + L +(x − x0 )m + n2!(m + n)! (2) = a 0 + a1(x − x0 ) + a 2(x − x0 )2 + L + a m + n(x − x0 )m + nĐể đơn giản ta coi x0 = 0. Ta cần tính các hệ số của Dn(x) và Qm(x) sao cho: Q (x) Tm + n (x) − m= 0 hay Tm+n(x)Dn(n) ‐ Qm(x) = 0 Dn (x)nghĩa là: (a 0 + a1x + L + a m + n x m + n )(1 + d1x + L + d n x n ) = (q 0 + q1x + L + q m x m ) (3) Cân bằng các số hạng cùng bậc ở hai vế ta có: ⎧a 0 = q 0⎪a + a d = q0 11⎪⎪ 1 ⎨a 2 + a1d1 + a 0d 2 = q 2⎪L⎪⎪⎩a m + a m −1d1 + a m −2d 2 + L + a m −nd n = q m ⎧a m +1 + a md1 + a m −1d 2 + L + a m −n+1d n = 0⎪a⎪ m + 2 + a m +1d1 + a md 2 + L + a m −n+ 2d n = 0 ⎨L⎪⎪⎩a m + n + a m + n −1d1 + a m + n−2d 2 + L + a md n = 0Trước hết ta giải (5) để tìm di và sau đó thay vào (4) để tìm qi. Ta xây dựng hàm padeapp() để tính xấp xỉ: function [num, den] = padeapp(f, xo, M, N, x0, xf) %Vao : f = Ham can xap xi trong doan [xo, xf] %Ra: num = Cac he so cua tu so % den = Cac he so cua mau so a(1) = feval(f, xo); h = .01; tmp = 1; for i = 1:M + N tmp = tmp*i*h; %i!h^i dix = difapx(i, [‐i i])*feval(f, xo + [‐i:i]*h)ʹ; %dao ham a(i + 1) = dix/tmp; %he so chuoi Taylor (4) (5) 225CuuDuongThanCong.com /> end for m = 1:N n = 1:N; A(m, n) = a(M + 1 + m ‐ n); b(m) = ‐a(M + 1 + m); end d = A\bʹ; %pt.(5) for m = 1: M + 1 mm = min(m ‐ 1,N); q(m) = a(m:‐1:m ‐ mm)*[1; d(1:mm)]; %pt.(4) end num = q(M + 1:‐1:1)/d(N); den = [d(N:‐1:1)ʹ 1]/d(N); %giam dan if nargout == 0 % ve ham thuc, khai trien taylor va ham Pade if nargin §7. NỘI SUY BẰNG ĐA THỨC HERMIT Trong một số trường hợp, ta cần tìm hàm đa thức khơng những đi qua các điểm cho trước mà cịn phải thoả mãn điều kiện về đạo hàm tại các điểm đó. Ta gọi đa thức như vậy là đa thức nội suy Hermit. Để đơn giản, ta khảo sát một đa thức bậc 3: h(x) = H 3 x 3 + H 2 x 2 + H1x + H0 (1) đi qua hai điểm (x0, y0), (x1, y1) và có các đạo hàm là y′0 , y′1 . Ta tìm các hệ số Hi bằng cách giải hệ phương trình: ⎧h(x 0 ) = H 3 x03 + H 2 x02 + H1x0 + H0 = y 0⎪32⎪h(x1 ) = H 3 x1 + H 2 x1 + H1x1 + H0 = y1 ⎨2′′h(x)=3Hx+2Hx+H=y03 02 010⎪⎪h′(x ) = 3H x 2 + 2H x + H = y′13 12 111⎩ (2) Các đạo hàm bậc nhất được tính gần đúng bằng: h(x0 + ε) − h(x0 ) y 2 − y 0=y′0 =εε (3) h(x1 ) − h(x1 − ε) y1 − y 3=y′1 =εεBây giờ ta tìm đa thưc nội suy Lagrange hay Newton đi qua 4 điểm: (x0, y0), (x 2 = x0 + ε ,y2 = y0 + y0′ ε) , (x 3 = x1 − ε , y3 = y1 − y′1ε) , (x1, y1) Hàm hermit() tạo nên phương trình (2): function H = hermit(x0, y0, dy0, x1, y1, dy1) A = [x0^3 x0^2 x0 1; x1^3 x1^2 x1 1; 3*x0^2 2*x0 1 0; 3*x1^2 2*x1 1 0]; b = [y0 y1 dy0 dy1]’; %Pt.(2) H = (A\b)’; Hàm hermits() dùng hàm hermit() để tính các hệ số của đa thức Hermit trên nhiều đoạn và giá trị nội suy: function [H,yi] = hermits(x, y, dy, xi) % Tim cac he so cua c da thuc Hermite tren c doan clc for n = 1:length(x)‐1 H(n,:) = hermit(0, y(n), dy(n), x(n + 1)‐x(n), y(n + 1), dy(n + 1)); 227CuuDuongThanCong.com /> end yi = ppval(mkpp(x, H),xi) Để nội suy ta dùng chương trình cthermite.m: clear all, clc x = [0 1 2 3]; y = [1 2 4 5]; dy = [0 2 4 6]; [h, y] = hermits(x, y, dy, 1.5) §8. BIẾN ĐỔI FOURIER 1. Biến đổi Fourrier: Tín hiệu thực tế thường bao gồm các thành phần có tần số khác nhau. Chuỗi Fourier và phép bíến đổi Fourier là cơng cụ tốn học dùng để phân tích đặc tính tần số của tín hiệu. Có 4 định nghĩa tương tự nhau về chuỗi và phép biến đổi Fourier, gồm: chuỗi Fourier liên tục theo t(CFS), phép biến đổi Fourier liên tục theo t(CFT), chuỗi Fourier gián đoạn theo t(DFS) và phép biến đổi Fourier gián đoạn theo t(DFT). Trong các cơng cụ này, DFT dễ dàng lập trình trên máy tính nên trong phần này ta sẽ chú ý đến nó. Giả sử chuỗi số liệu { x[n] = x(nT), n = 0 : M ‐ 1} với T là chu kì lấy mẫu có được bằng cách lấy mẫu một tín hiệu liên tục x(t) T lần trong một giây. N cặp điểm DFT và iDFT được định nghĩa bằng: N −1 DFT: X(k) = ∑ x[n]e ‐j2πnk/N (1a) n =0N −11X(k)e j2πnk/N (1b) ∑N n =0Nói chung hệ số DFT của X(k) là một số phức và nó xác định biên độ và pha của thành phần tín hiệu có tần số số Ωk = kΩ0(rad), tương ứng với tần số tương tự ωk = kω0 = kΩ0/T = 2πk/NT (rad/s). Ta gọi Ω0 = 2π/N và ω0 = 2π/NT là các tần số cơ bản số và tương tự (tần số phân giải) vì đây là hiệu tần số có thể phân biệt bởi N điểm DFT. DFT và DFS có cùng bản chất nhưng khác nhau về phạm vi thời gian/tần số. Cụ thể là tín hiệu x[n] và DFT X[k] của nó kéo dài hữu hạn trên phạm vi thời gian/tần số {0 ≤ n ≤ N‐1} và {0 ≤ k ≤ N‐1}. Tín hiệu x[n] được iDFT: x[n] =228CuuDuongThanCong.com /> phân tích bởi DFS và DFS của nó X(k) là chu kì tín hiệu với chu kì N trên tồn bộ tập số ngun. Biến đổi Fourier nhanh FFT là thuật tốn hiệu quả để tính DFT và iDFT được xây dựng bằng cách dùng tính chu kì và tính đối xứng cuả nhân tử ei2πnk/N để giảm bớt số nhân tử phức từ N2 thành (N/2)log2N )N thể hiện kích thước của DFT. Hàm MATLAB fft() và ifft() thực hiện thuật tốn đối với N = 2l (l là số ngun khơng âm). Nếu độ dài M của chuỗi số liệu ban đầu khơng phải là bội số của 2, có thể mở rộng bằng cách đệm thêm số 0 vào cuối chuỗi và gọi là đệm zero. Ta xem xét hiệu qủa này bằng cách thực hiện đoạn lệnh trong ctcompdftfft.m. %So sanh phep bien doi Fourier nhanh va roi rac clear, clf N = 2^10; n = [0:N ‐ 1]; x = cos(2*pi*200/N*n)+ 0.5*sin(2*pi*300/N*n); tic %ngung dong ho for k = 0:N ‐ 1 X(k+1) = x*exp(‐j*2*pi*k*n/N).ʹ; end %DFT k = [0:N ‐ 1]; for n = 0:N ‐ 1 xr(n + 1) = X*exp(j*2*pi*k*n/N).ʹ; end %IDFT time_dft = toc plot(k,abs(X)) pause, hold on tic X1 = fft(x); %FFT xr1 = ifft(X1); %IFFT time_fft = toc %dua ra thoi gian thuc hien clf, plot(k,abs(X1),ʹrʹ) %pho bien do Chạy đoạn lệnh và so sánh thời gian thực hiện 1024 điểm tính DFT/iDFT và FFT/iFFT. 229CuuDuongThanCong.com /> 2. Ý nghĩa vật lý của biến đổi Fourrier rời rạc: Để hiểu được ý nghĩa vật lí của FFt ta thực hiện các lệnh trong chương trình ctmeanning.m. Chương trình cho ta phổ biên độ của tín hiệu x(t) = sin(1.5πt) + 0.5cos(3πt) (2) được lấy mẫu mỗi T s. Từ các kết quả ta thấy khi T = 0.1 và N = 32 thì Xa(k) lớn tại k = 2 và k= 5. Lúc đó kω0 = 2πk/NT = 2πk/3.2 ≈ 1.5π và 3.125π ≈ 3π. Khi T = 0.05 và N = 64 thì Xb(k) cũng lớn tại k = 2 và k = 5. Lúc đó kω0 = 1.25π ≈ 1.5π và 3.125π ≈ 3π. Khi T = 0.1 và N = 64 thì Xc(k) lớn tại k = 4 ,k = 5, k = 9 và k = 10. Lúc đó kω0 = 2πk/NT = 2πk/6.4 ≈ 1.25π ~ 1.5625π và 2.8125π ~ 3π. Khi T = 0.1 và N = 64 thì Xd(k) lớn tại k = 5 và k = 10. Lúc đó kω0 = 1.5625π ≈ 1.5π và 3.125π ≈ 3π. Tồn tại nhiều phổ DFT khác nhau của cùng một tín hiệu tương tự, tuỳ thuộc vào kích thước DFT, chu kì lấy mẫu, khoảng biến thiên của hàm và đệm zero. So sánh với phổ tại T = 0.1s, phổ tại T = 0.05s có phạm vi tần số tương tự [0, 2π/Tb] rộng hơn nhưng có cùng tần số phân giải tương tự là ω0 = Ω0/Tb = 2π/NbTb = π/1.6 = 2π/NaTa. Phổ khi có đệm zero trơn. clear, clf w1 = 1.5*pi; w2 = 3*pi; N = 32; n = [0:N ‐ 1]; T = 0.1; %chu ki lay mau t = n*T; xan = sin(w1*t) + 0.5*sin(w2*t); subplot(421) stem(t,xan,ʹ.ʹ) k = 0:N ‐ 1; Xa = fft(xan); dscrp=norm(xan‐real(ifft(Xa))) subplot(423) stem(k,abs(Xa),ʹ.ʹ) N = 32; n = [0:N ‐ 1]; 230CuuDuongThanCong.com /> T = 0.1; t = n*T; xan = sin(w1*t) + 0.5*sin(w2*t); subplot(422) stem(t,xan,ʹ.ʹ) k = 0:N ‐ 1; Xa = fft(xan); Dscrp = norm(xan ‐ real(ifft(Xa))) subplot(424) stem(k, abs(Xa),ʹ.ʹ) N = 64; n = [0:N ‐ 1]; T = 0.05; t = n*T; xbn = sin(w1*t) + 0.5*sin(w2*t); subplot(425) stem(t,xbn,ʹ.ʹ) k = 0:N ‐ 1; Xb = fft(xbn); subplot(427) stem(k,abs(Xb),ʹ.ʹ) N = 64; n = [0:N‐1]; T = 0.1; t = n*T; xbn = sin(w1*t) + 0.5*sin(w2*t); subplot(426) stem(t, xbn,ʹ.ʹ) k = 0:N ‐ 1; Xb = fft(xbn); subplot(428) stem(k, abs(Xb),ʹ.ʹ) Ta có nhiều phổ DFT cho cùng một tín hiệu tương tự, tuỳ thuộc vào kích thước DFT, chu kì lấy mẫu, khoảng lấy mẫu và đệm zero. So sánh phố khi giảm chu kì lấy mẫu T từ 0.1s đến 0.05s 231CuuDuongThanCong.com /> 3. Nội suy bằng các dùng biến đổi Fourrier rời rạc: Ta dùng DFS/DFT để nội suy dãy x[n] nhận được từ kết quả lấy mẫu tín hiệu ở khoảng cách cách đều. Thủ tục gồm hai bước: lấy N điểm FFT X(k) của x[n] và dùng cơng thức: 1j2 πkt / NT%xˆ (t) =X(k)e∑N |k| (5) N / 2 −11⎧⎫= ⎨X(0) + 2 ∑ Real ⎣⎡ X(k)e j2 πkt / NT ⎦⎤ + X(N / 2)cos(π/T) ⎬N⎩k =1⎭Ta xây dựng hàm nội suy interpdfs(): function [xi,Xi] = interpdfs(T, x, Ws, ti) %T : chu li lay mau %x : thu tu roi rac hoa %Ws: tan so dung chuan (1.0 = pi[rad]) %ti: khoang thoi gian noi suy if nargin w1 = pi; w2 = .5*pi; %hai tan so N = 32; n = [0:N ‐ 1]; T = 0.1; t = n*T; x = sin(w1*t)+0.5*sin(w2*t)+(rand(1,N) ‐ 0.5); %0.2*sin(20*t); ti = [0:T/5:(N ‐ 1)*T]; subplot(411), plot(t,x,ʹk.ʹ) %so lieu ban dau title(ʹSo lieu ban dau va ket qua noi suyʹ) [xi,Xi] = interpdfs(T,x,1,ti); hold on, plot(ti,xi,ʹrʹ) %tai tao tin hieu k = [0:N ‐ 1]; subplot(412), stem(k,abs(Xi),ʹk.ʹ) %pho ban dau title(ʹPho ban dauʹ) [xi,Xi] = interpdfs(T,x,1/2,ti); subplot(413), stem(k,abs(Xi),ʹr.ʹ) %pho da loc title(ʹPho da locʹ) subplot(414), plot(t,x,ʹk.ʹ, ti,xi,ʹrʹ) %tin hieu da loc title(ʹTin hieu da locʹ) §9. XẤP XỈ HÀM BẰNG PHƯƠNG PHÁP BÌNH PHƯƠNG BÉ NHẤT 1. Khái niệm chung: Trong các mục trước ta đã nội suy giá trị của hàm. Bài tốn đó là cho một hàm dưới dạng bảng số và phải tìm giá trị của hàm tại một giá trị của đối số khơng nằm trong bảng. Trong thực tế, bên cạnh bài tốn nội suy ta cịn gặp một dạng bài tốn khác. Đó là tìm cơng thức thực nghiệm của một hàm. Nội dung bài tốn là từ một loạt các điểm cho trước (có thể là các giá trị của một phép đo nào đó) ta phải tìm một hàm xấp xỉ các giá trị đã cho. Ta sẽ dùng phương pháp bình phương tối thiểu để giải bài tốn. Giả sử có mẫu quan sát (xi, yi) của hàm y = f(x). Ta chọn hàm f(x) có dạng: f(x) = a0f0(x) + a1f1(x) + a2f2(x)... (1) Trong đó các hàm f0(x), f1(x), f2(x) v.v. là (m+1) hàm độc lập tuyến tính mà ta có thể chọn tuỳ ý và các hệ số ai là tham số chưa biết mà ta phải xác định dựa 233CuuDuongThanCong.com /> vào hệ hàm đã chọn và các điểm quan sát. Sai số giữa trị đo được và trị tính theo (1) là : ei = yi ‐ f(xi) (2) Sai số này có thể âm hay dương tuỳ từng giá trị của yi. Khi dùng phương pháp bình phương bé nhất ta xét bình phương của sai số tại một điểm: e i2 = [ y i − f(x i )] 2 (3) Với n điểm tổng bình phương của sai số sẽ là : nni =1i =1S = ∑ e i2 = ∑ {y i − [a 0 f0 (x i ) + a1f1 (x i ) + ⋅ ⋅ ⋅ + a m fm (x i )]} 2Rõ ràng S là hàm của các giá trị cần tìm ai và chúng ta sẽ chọn các ai sao ∂Sphải bằng khơng. cho S đạt giá trị min, nghĩa là các đạo hàm ∂a iTa sẽ xét các trường hợp cụ thể. 2. Hàm xấp xỉ có dạng đa thức: Trong trường hợp tổng qt ta chọn hệ hàm xấp xỉ là một đa thức, nghĩa là: f(x) = a0 + a1x + a2x2 +∙∙∙+ amxm Vậy hàm S là : S = ( y i − a 0 + a1x + a 2 x 2 + ⋅⋅⋅ + a m x m ) Theo điều kiện đạo hàm 2 ∂S= 0 ta nhận được hệ phương trình: ∂a inn⎧a n x m + am −1+⋅⋅⋅+=xnayi∑∑∑m −1i0⎪ m i =1 ii =1i =1⎪ nnnn⎪a ∑ x m + 1 + a ∑ x m + ⋅ ⋅ ⋅ +a ∑ x = ∑ x ym −1i0ii i⎪ m i =1 ii =1i =1i =1⎪ nnnnm+2m +122⎪⎪a m ∑ xi + a m −1 ∑ x i + ⋅ ⋅ ⋅ +a 0 ∑ xi = ∑ x i y i i =1i =1i =1⎨ i =1nnn⎪ n m+3m+23++⋅⋅⋅+=axaxax∑ xi3 y i⎪ m∑ im −1 ∑ i0∑ ii =1i =1i =1⎪ i =1⎪⋅ ⋅ ⋅⎪ nnnn⎪a m ∑ xi2 m + a m −1 ∑ x i2 m −1 + ⋅ ⋅ ⋅ +a 0 ∑ x im = ∑ x im y i⎪⎩ i =1i =1i =1i =1 Đây là một hệ phương trình tuyến tính. Giải nó ta nhận được các gía trị ai. Ta xây dựng hàm polynomfit() thực hiện thuật tốn trên: 234CuuDuongThanCong.com />

Tài liệu liên quan

Kế toán tài chính - chương 4
- 39
- 984
- 1
Tài liệu Tài chính tiền tệ_ Chương 4 tiền tệ và lưu thông tiền tệ pdf
- 2
- 598
- 1
Tài liệu CHƯƠNG 3: NỘI SUY VÀ XẤP XỈ HÀM ppt
- 31
- 807
- 5
quản trị hành chính văn phòng - chương 4-nâng cao hiệu quả công việc hành chính văn phòng
- 20
- 1
- 1
Chương 5 NỘI SUY VÀ XẤP XỈ HÀM pptx
- 7
- 3
- 36
Chương 5 - Nội suy và xấp xỉ hàm doc
- 24
- 606
- 2
Tài liệu quản lý vận hành - Chương 4 Dự báo
- 124
- 948
- 3
chương 4 phát triển và các vấn đề về tài nguyên thiên nhiên
- 24
- 811
- 1
Chương 3: Nội suy và xấp xỉ hàm số doc
- 34
- 1
- 3
Kinh tế tài chính - Chương 4 - Kế toán TSCĐ, BĐS đầu tư và đầu tư XDCB pot
- 75
- 613
- 0