Giáo trình Turbo C nâng cao và C++ - Chương 11: Nội suy và xấp xỉ hàm

180 Ch−ơng 11 : nội suy và xấp xỉ hàm Đ1.Nội suy Lagrange Trong thực tế nhiều khi phải phục hồi một hàm y = f(x) tại mọi 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 toán.Vì vậy nảy sinh vấn đề toá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.Bây giờ ta cần tìm đa thức : Pn(x) = aox n + a1x n-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,luôn có đạo hàm và nguyên hàm.Việc tính giá trị của nó theo thuật toá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 : )xx)...(xx)(xx)...(xx( )xx)...(xx)(xx)...(xx( L ni1ii1ii0i n1i1i0 i −−−− −−−−= +− +− Rõ ràng là Li(x) là một đa thức bậc n và : ⎪⎩ ⎪⎨ ⎧ = ≠ = ij0 ij1 )x(L ji 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 0i iin )x(L)x(f)x(P Ta 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 y1 Đa thức nội suy sẽ là : P1(x) = yoL0(x) + y1L1(x1) 10 1 0 xx xx L − −= 01 0 1 xx xx L − −= nên 01 0 1 10 1 01 xx xx y xx xx y)x(P − −+− −= Nh− 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 y0 y1 y2 Đa thức nội suy sẽ là : P2(x) = yoL0(x) + y1L1(x1) + y2L2(x2) )xx)(xx( )xx)(xx( L 2010 21 0 −− −−= )xx)(xx( )xx)(xx( L 2101 20 1 −− −−= 181 )xx)(xx( )xx)(xx( L 1202 10 2 −− −−= Nh− vậy P1(x) là một đa thức bậc hai đối với x Trên cơ sở thuật toán trên ta có ch−ơng trình tìm đa thức nội suy của một hàm khi cho tr−ớc các điểm và sau đó tính trị số của nó tại một giá trị nào đó nh− sau : Ch−ơng trình 11-1 #include #include #include #define max 21 int maxkq,n; float x[max],y[max],a[max],xx[max],yy[max]; float x0,p0; void main() { int i,k; char ok ; void vaosolieu(void); float lagrange(int,float [],float [],float); void inkq(void); clrscr(); printf("%24cNOI SUY DA THUC LAGRANGE\n",' '); vaosolieu(); k=0; ok='c'; while (ok=='c') { printf("Tinh gia tri cua y voi x la x0 = "); scanf("%f",&x0); p0=lagrange(n,x,y,x0); printf("Gia tri cua y = %15.5f\n",p0); printf("\n"); k=k+1; maxkq=k; xx[k]=x0; yy[k]=p0; flushall(); printf("Tinh tiep khong(c/k)?"); scanf("%c",&ok); } inkq(); 182 } void vaosolieu() { int i,t; char ok; printf("\n"); printf("Ham y = f(x)\n"); printf("So cap (x,y) nhieu nhat la max = 20\n"); printf("So diem da cho truoc n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } printf("\n"); printf(" SO LIEU BAN VUA NHAP\n"); printf(" x y\n"); for (i=1;i<=n;i++) printf("%8.4f %8.4f\n",x[i],y[i]); ok=' '; t=1; flushall(); while (t) { printf("\nCo sua so lieu khong(c/k):?"); scanf("%c",&ok); if (toupper(ok)=='C') { printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); flushall(); } if (toupper(ok)!='C') t=0; } } float lagrange(int n,float x[max],float y[max],float x0) { int i,k; 183 float g0; p0=0.0; for (k=1;k<=n;k++) { g0=1.0; for (i=1;i<=n;i++) if (i!=k) g0=g0*(x0-x[i])/(x[k]-x[i]); p0=p0+y[k]*g0; } return(p0); } void inkq() { int i,j,k; printf("\n"); printf("%24cBANG SO LIEU\n",' '); printf("%18cx %24cy\n",' ',' '); for (i=1;i<=n;i++) printf("%20.4f %25.4f\n",x[i],y[i]); printf("\n"); printf("%24cKET QUA TINH TOAN\n",' '); printf("%14cx %10cy\n",' ',' '); for (k=1;k<=maxkq;k++) printf("%15.5f %15.5f\n",xx[k],yy[k]); getch(); } Giả sử ta có bảng các giá trị x,y : x 0 3 -2 2 4 y 0 -3.75 10 -2 4 vậy theo ch−ơng trình tại x = 2.5 y = -3.3549. Đ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 ts đ−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à : 184 ji ji ji xx yy ]x,x[y − −= Tỉ hiệu cấp hai của y tại xi,xj,xk là : ki kjji kji xx ]x,x[y]x,x[y ]x,x,x[y − −= v.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 : 0 0nn 0n xx )x(P)x(P ]x,x[P − −= là một đa thức bậc (n-1).Tỉ hiệu cấp 2 tại x,x0,x1 : 1 10n0n 10n xx ]x,x[P]x,x[P ]x,x,x[P − −= là một đa thức bậc (n-2) v.v và tới tỉ hiệu cấp (n+1) thì : 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).Ngoà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ó : h y ]x,x[y 010 ∆= h2 y ]x,x,x[y 20 2 210 ∆= ........... 185 )h!n( y ]x.,.,.x[y n0 n n0 ∆= Bây giờ đặt x = x0 + ht trong đa thức Newton tiến ta đ−ợc : y !n )1nt.(.).1t(t y !2 )1t(t yty)htx(P 0 n 0 2 000n ... ∆+−−∆−∆ ++++=+ 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ó : y 2 )1t(t yty)htx(P 0 2 0002 ∆−∆ ++=+ 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 : y !n )1nt.(.).1t(t y !2 )1t(t yty)htx(P n n n 2 nn0n ... ∇−++∇+∇ ++++=+ Ví dụ : Cho hàm nh− bảng sau : x 0.1 0.2 0.3 0.4 y 0.09983 0.19867 0.29552 0.38942 Ta tính giá trị của hàm tại 0.14 bằng đa thức nội suy Newton vì các mốc cách đều h = 0.1.Ta có bảng sai phân sau : i x y ∆y ∆2y ∆3y 0 0.1 0.09983 0.09884 1 0.2 0.19867 - 0.00199 0.09685 -0.00096 2 0.3 0.29552 - 0.00295 0.09390 3 0.4 0.38942 Ta dùng công thức Newton tiến với điểm gốc là x0 = 0.1.h = 0.1.Với x = 0.14 ta có 0.14 = 0.1 + 0.1t nên t = 0.4 và kết quả là : 13954336.000096.0 !3 )2t)(1t(t 00199.0 !2 )1t(t 099884.0.t09983.0)t1.01.0(P =−−−−++=+ Ch−ơng trình nội suy Newton nh− sau : Ch−ơng trình 11-2 186 //Noi suy Newton #include #include #include #define max 11 void main() { int i,j,k,n,t; float a[max],b[max],x[max],y[max]; char ok; float x0,p; clrscr(); printf("So diem da cho n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } printf("%10cBANG SO LIEU\n",' '); printf("%8cx%30cy\n",' ',' '); for (i=1;i<=n;i++) printf("%4c%8.4f%23c%8.4f\n",' ',x[i],' ',y[i]); ok=' '; t=0; flushall(); while (t) { printf("Co sua so lieu khong(c/k): "); scanf("%c",&ok); if (toupper(ok)=='C') { printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); flushall(); } if (toupper(ok)!='C') t=0; } a[1]=y[1]; for (j=1;j<=n-1;j++) 187 { for (i=1;i<=n-j;i++) y[i]=(y[i+1]-y[i])/(x[i+j]-x[i]); a[j+1]=y[1]; } b[n]=a[n]; for (k=n-1;k>=1;k--) { for (j=n-1;j>=1;j--) b[j]=a[j] ; for (i=n-1;i>=k;i--) a[i]=a[i]-b[i+1]*x[k]; } for (i=n;i>=1;i--) printf("He so bac %d la :%8.4f\n",i-1,a[i]); printf("\n"); k=0; ok='c'; flushall(); while (ok=='c') { printf("Tinh gia tri cua y tai x = "); scanf("%f",&x0); p=0; for (k=n;k>=1;k--) p=p*x0+a[k]; printf("Tri so noi suy tai x0 = %4.2f la : %10.5f\n",x0,p); getch(); printf("Ban co muon tinh tiep cac diem khac khong(c/k)"); do scanf("%c",&ok); while ((ok!='c')&&(ok!='k')); } } Dùng ch−ơng trình này nội suy các giá trị cho trong bảng sau 0 0.2 0.4 0.6 0.8 1.0 1 1.2214027 6 1.4918247 1.8221188 2.2255409 3 2.7182818 3 ta có các hệ số của đa thức nội suy : 0.0139(bậc 5),0.0349(bậc 4),0.1704(bậc3),0.4991(bậc 2),1.0001(bậc 1) và 1.0000(bậc 0). Đ3.Nội suy Aitken Một dạng khác của đa thức nội suy đ−ợc xác định bằng thuật toán Aitken.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 : 188 01 11 00 01 xx xxy xxy )x(P − − − = là một đa thức bậc 1 : 01 0 1 10 1 001 xx xx y xx xx y)x(P − −+− −= .Khi x = x0 thì : 0 01 011 000 001 yxx xxy xxy )x(P =− − − = Khi x = x1 thì : 1 01 111 100 101 yxx xxy xxy )x(P =− − − = Đa thức nội suy Lagrange của f(x) qua 3 điểm x0,x1,x2 có dạng : 02 212 001 012 xx xx)x(P xx)x(P )x(P − − − = và là một đa thức bậc 2: )xx)(xx( )xx)(xx( y )xx)(xx( )xx)(xx( y )xx)(xx( )xx)(xx( y)x(P 1202 10 2 2101 20 1 2010 21 0012 −− −−+−− −−+−− −−= Khi x = x0 thì : 0 02 0212 000 0012 yxx xx)x(P xxy )x(P =− − − = Khi x = x1 thì : 1 02 121 101 1012 yxx xxy xxy )x(P =− − − = Khi x = x2 thì : 2 02 222 20201 2012 yxx xxy xx)x(P )x(P =− − − = Tổng quát đa thức nội suy Lagrange qua n điểm là : 02 nn..12 0)1n..(01 n..012 xx xx)x(P xx)x(P )x(P − − − = − Nh− 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 toán nh− vậy gọi là sơ đồ Neville-Aitken. Ví dụ : Cho các cặp điểm (0,0.4),(1.4,1.5),(2.6,1.8),(3.9,2.6),tính y tại x=2 189 97143.1 04.1 6.05.1 24.0 xx xxy xxy )x(P 01 11 00 01 =− − − =− − − = 65.1 4.16.2 6.08.1 6.05.1 xx xxy xxy )x(P 12 22 11 12 =− − =− − − = 7242.1 06.2 6.065.1 297143.1 xx xx)x(P xx)x(P )x(P 02 212 001 012 =− − =− − − = 4308.1 6.29.3 9.16.2 6.08.1 xx xxy xxy )x(P 23 33 22 23 =−=− − − = 5974.1 4.19.3 9.14308.1 6.065.1 xx xx)x(P xx)x(P )x(P 13 323 112 123 =− − =− − − = 6592.1 09.3 9.15974.1 27242.1 xx xx)x(P xx)x(P )x(P 03 3123 0012 0123 =− − =− − − = Ch−ơng trình đ−ợc viết nh− sau Ch−ơng trình 11-3 //Noi suy Aitken #include #include #include #define max 11 void main() { float x[max],y[max],yd[max]; float x1; int j,k,n,n1; clrscr(); printf("Cho so diem da co n = "); scanf("%d",&n1); n=n1-1 ; for (k=0;k<=n;k++) { printf("x[%d] = ",k+1); scanf("%f",&x[k]); printf("y[%d] = ",k+1); scanf("%f",&y[k]); } printf("Cho diem can tinh gia tri cua ham x1 = "); 190 scanf("%f",&x1); for (k=0;k<=n-1;k++) { yd[k]=(y[k]*(x1-x[k+1])-y[k+1]*(x1-x[k]))/(x[k]-x[k+1]); if (k!=0) for (j=k-1;j>=0;j--) yd[j]=(yd[j]*(x1-x[k+1])-yd[j+1]*(x1-x[j]))/(x[j]-x[k+1]); } printf("Gia tri ham tai x = %6.3f la y = %8.4f\n",x1,yd[0]); getch(); } Dùng ch−ơng trình này để nội suy các cặp số (1,3),(2,5),(3,7),(4,9) và (5,11) tại x = 2.5 ta có y = 6. Đ4.Xấp xỉ hàm bằng ph−ơng pháp bình ph−ơng bé nhất Trong các mục tr−ớc ta đã nội suy giá trị của hàm.Bài toá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 toán nội suy ta còn gặp một dạng bài toá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 toá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 toá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 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 : [ ]2ii2i )x(fye −= (3) Với n điểm tổng bình ph−ơng của sai số sẽ là : [ ]{ }∑∑ == +⋅⋅⋅++−== n 1i 2 inni11i00i n 1i 2 i )x(fa)x(fa)x(fayeS Rõ 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 cho S đạt giá trị min,nghĩa là các đạo hàm ia S ∂ ∂ phải bằng không.Ta sẽ xét các tr−ờng hợp cụ thể. 1.Hàm xấp xỉ có dạng đa thức : Trong tr−ờng hợp tổng quát ta chọn hệ hàm xấp xỉ là một đa thức,nghĩa là : f(x) = a0 + a1x + a2x 2 +...+ amx m Vậy hàm S là : ( )2m210i xaxaxaayS +⋅⋅⋅+++−= Theo điều kiện đạo hàm 0 a S i =∂ ∂ ta nhận đ−ợc hệ ph−ơng trình: 191 ⎪⎪ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎪⎪ ⎨ ⎧ =+⋅⋅⋅++ ⋅⋅⋅ =+⋅⋅⋅++ =+⋅⋅⋅++ =+⋅⋅⋅++ =+⋅⋅⋅++ ∑∑∑ ∑ ∑∑∑ ∑ ∑∑∑ ∑ ∑∑∑ ∑ ∑∑ ∑ === = − − === = + − + === = + − + === = − + == = − − n 1i i m i n 1i m i0 n 1i n 1i 1m2 i1m m2 im n 1i i 3 i n 1i 3 i0 n 1i n 1i 2m i1m 3m im n 1i i 2 i n 1i 2 i0 n 1i n 1i 1m i1m 2m im n 1i ii n 1i i0 n 1i n 1i m i1m 1m im n 1i i0 n 1i n 1i 1m i1m m im yxxaxaxa yxxaxaxa yxxaxaxa yxxaxaxa ynaxaxa Đâ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.Sau đây là ch−ơng trình viết theo thuật toán trên. Ch−ơng trình 11-4 //Xap xi da thuc #include #include #include #define max 11 void main() { int i,j,k,m,n,p,kp,t; float a[max],x[max],y[max],y1[max]; float b[max][max]; char ok; float s,sx,s1,c,d; clrscr(); printf("PHUONG PHAP BINH PHUONG TOI THIEU"); printf("\n"); printf("Cho bac cua da thuc xap xi m = "); scanf("%d",&m); printf("So diem da cho n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } x[0]=1; printf("\n"); printf("%4cBANG SO LIEU\n",' '); 192 printf("%8cx%30cy\n",' ',' '); for (i=1;i<=n;i++) printf("%4c%8.4f%20c%8.4f\n",' ',x[i],' ',y[i]); ok=' '; t=1; flushall(); while (t) { printf("Co sua so lieu khong(c/k): "); scanf("%c",&ok); if (toupper(ok)=='C') { printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); flushall(); } if (toupper(ok)!='C') t=0; } //for (i=0;i<=n;i++) //a[i]=0.0; printf("\n"); printf("CAC GIA TRI DA CHO"); printf("\n"); printf("X = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',x[i]); printf("\n"); printf("Y = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',y[i]); printf("\n"); for (p=0;p<=m;p++) { y1[p]=0.0; for (i=1;i<=n;i++) { sx=1.0; for (j=1;j<=p;j++) sx*=x[i]; y1[p]+=y[i]*sx; } } for (p=0;p<=m;p++) for (k=0;k<=m;k++) 193 { kp=k+p; b[p][k]=0.0; for (i=1;i<=n;i++) { sx=1.0; for (j=1;j<=kp;j++) sx*=x[i]; b[p][k]+=sx; } } for (i=0;i<=m-1;i++) { c=1.0/b[i][i]; for (k=i+1;k<=m;k++) { d=b[i][k]; for (j=i+1;j<=m;j++) b[k][j]-=b[i][j]*c*d; y1[k]-=y1[i]*c*d; b[i][k]*=c; } y1[i]*=c; } y1[m]/=b[m][m]; for (i=m-1;i>=0;i--) for (j=i+1;j<=m;j++) y1[i]-=b[i][j]*y1[j]; printf("\n"); printf("CAC HE SO CUA DA THUC CAN TIM"); printf("\n"); for (i=0;i<=m;i++) printf("a[%d] = %10.5f\n",i,y1[i]); getch(); } Với các giá trị x,y đo đ−ợc theo bảng x 7 8 9 10 11 12 13 y 7,4 8,4 9,1 9,4 9,5 9,5 9,4 ta có n = 7 và chọn m = 2 và tính đ−ợc theo ch−ơng trình các hệ số : a0 = -0.111905 ; a1 = 2.545238 ; a2 = -4.857143 và hàm xấp xỉ sẽ là : f(x) = -0.111905 + 2.545238x -4.857143x2 2.Hàm dạng Aecx : Khi các số liệu thể hiện một sự biến đổi đơn điệu ta dùng hàm xấp xỉ là y = Aecx.Lấy logarit hai vế ta có : lny = lnA + cxlne 194 Theo điều kiện đạo hàm 0 a S i =∂ ∂ ta có hệ ph−ơng trình : ⎪⎪⎩ ⎪⎪⎨ ⎧ =+ =+ ∑ ∑∑ ∑ ∑ = == = = n 1i n 1i ii n 1i i 2 i n 1i n 1i ii ylnxxAlnxc ylnAlnnxc Giải hệ ph−ơng trình này ta có các hệ số A và c : Ch−ơng trình 11-5 //xap_xi_e_mu; #include #include #include #include #define max 11 void main() { int i,n,t; float x[max],y[max]; char ok; float a,b,c,d,e,f,d1,d2,d3; clrscr(); printf("PHUONG PHAP BINH PHUONG TOI THIEU"); printf("\n"); printf("So diem da cho n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } x[0]=1.0; printf("%4cBANG SO LIEU\n",' '); printf("%8cx%30cy\n",' ',' '); for (i=1;i<=n;i++) printf("%4c%8.4f%23c%8.4f\n",' ',x[i],' ',y[i]); ok=' '; t=1; while (t) { printf("Co sua so lieu khong(c/k): "); scanf("%c",&ok); if (toupper(ok)=='C') { 195 printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } if (toupper(ok)!='C') t=0; } printf("CAC GIA TRI DA CHO"); printf("\n"); printf("X = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',x[i]); printf("\n"); printf("Y = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',y[i]); printf("\n"); a=0.0; for (i=1;i<=n;i++) a+=x[i]; b=n; c=0.0; for (i=1;i<=n;i++) c+=log(y[i]); d=0.0; for (i=1;i<=n;i++) d+=x[i]*x[i]; e=0.0; for (i=1;i<=n;i++) e+=x[i]*log(y[i]); d1=a*a-d*b; d2=c*a-e*b; d3=a*e-c*d; c=d2/d1; a=d3/d1; printf("\n"); printf("He so A = %8.4f",exp(a)); printf(" va so mu c = %8.4",c); printf("\n"); printf("\nBANG CAC GIA TRI TINH TOAN"); printf("\n"); printf("%5cx%28cy\n",' ',' '); for (i=1;i<=n;i++) { printf("%8.4f%21c%8.4f\n",x[i],' ',exp(a)*exp(c*x[i])); } 196 getch(); } Với các giá trị x,y đo đ−ợc theo bảng x 0 2 4 6 8 10 12 y 128 0 635 324 162 76 43 19 ta có n = 7 và tính đ−ợc theo ch−ơng trình các hệ số : A = 1285.44 va c = -0.3476 và hàm xấp xỉ sẽ là : f(x) = 1285.44 3.Hàm dạng Axq : Khi các số liệu thể hiện một sự biến đổi đơn điệu ta cũng có thể dùng hàm xấp xỉ là y = Axq.Lấy logarit hai vế ta có : lny = lnA + qlnx Theo điều kiện đạo hàm triệt tiêu ta có hệ ph−ơng trình : ⎪⎪⎩ ⎪⎪⎨ ⎧ =+ =+ ∑ ∑∑ ∑ ∑ = == = = n 1i n 1i ii n 1i ii 2 n 1i n 1i ii ylnxlnxlnAlnxlnq ylnAlnnxlnq Giải hệ ph−ơng trình này ta có các hệ số A và q : Ch−ơng trình 11-6 //xap_xi_x_mu; #include #include #include #include #define max 11 void main() { int i,n,t; float x[max],y[max]; char ok; float a,b,c,d,e,f,d1,d2,d3; clrscr(); printf("PHUONG PHAP BINH PHUONG TOI THIEU"); printf("\n"); printf("So diem da cho n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); 197 } x[0]=1.0; printf("%4cBANG SO LIEU\n",' '); printf("%8cx%30cy\n",' ',' '); for (i=1;i<=n;i++) printf("%4c%8.4f%23c%8.4f\n",' ',x[i],' ',y[i]); ok=' '; flushall(); t=1; while (t) { printf("Co sua so lieu khong(c/k): "); scanf("%c",&ok); if (toupper(ok)=='C') { printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[",i,"] = "); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } if (toupper(ok)!='C') t=0; } printf("\n"); printf("\nCAC GIA TRI DA CHO"); printf("\n"); printf("X = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',x[i]); printf("\n"); printf("Y = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',y[i]); printf("\n"); a=0.0; for (i=1;i<=n;i++) a+=log(x[i]); b=n; c=0.0; for (i=1;i<=n;i++) c+=log(y[i]); d=0.0; for (i=1;i<=n;i++) d+=log(x[i])*log(x[i]); e=0.; for (i=1;i<=n;i++) e+=log(x[i])*log(y[i]); 198 d1=a*a-d*b; d2=c*a-e*b; d3=a*e-c*d; c=d2/d1; a=d3/d1; printf("\n"); printf("He so A = %8.4f",exp(a)); printf(" va so mu q = %8.4f\n",c); printf("\n"); printf("\nBANG CAC GIA TRI TINH TOAN\n"); printf("%5cx%27cy\n",' ',' '); for (i=1;i<=n;i++) { printf("%8.4f%20c%8.4f\n",x[i],' ',exp(a)*exp(c*log(x[i]))); } getch(); } Với các giá trị x,y đo đ−ợc theo bảng x 1 2 4 5 6 y 7.1 27.8 62.1 110 161 ta có n = 5 và tính đ−ợc theo ch−ơng trình các hệ số : A = 7.1641 và q = 1.9531 và hàm xấp xỉ sẽ là : f(x) = 1285.44x1.9531 4.Hàm l−ợng giác : Khi quan hệ y=f(x) có dạng tuần hoàn ta dùng hàm xấp xỉ là tổ hợp tuyến tính của các hàm sin và cosin dạng : ∑ ∑ = = ω+ω+= n 1i n 1i ii0 )xisin(b)xicos(aa)x(f Để đơn giản tr−ớc hết ta xét hàm chỉ có một số hạng sin-cos,nghĩa là : xsinbxcosaa)x(f 110 ω+ω+= Hàm S sẽ có dạng : [ ]∑ = ω+ω+−= n 1i 2 110i )xsinbxcosaa(yS Theo điều kiện đạo hàm triệt tiêu ta có hệ ph−ơng trình đối với các hệ số dạng : ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ω ω=⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ωωωω ωωωω ωω ∑∑ ∑ ∑∑∑ ∑∑∑ ∑∑ xsiny xcosy y b a a xsinxsinxcosxsin xsinxcosxcosxcos xsinxcosn 1 1 0 2 2 Do : 0 n xsinxcos 2 1 n xcos 2 1 n xsin 0 n xcos 0 n xsin 22 =ωω =ω=ω =ω=ω ∑ ∑∑ ∑∑ nên hệ ph−ơng trình có dạng đơn giản : 199 ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ ω ω=⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ⎥⎥⎦ ⎤ ⎢⎢⎣ ⎡ ∑∑ ∑ xsiny xcosy y b a a 2n00 02n0 00n 1 1 0 Giải hệ ta có : ∑∑∑ ω=ω== xsinyn2bxcosyn2an y a 110 Trong tr−ờng hợp tổng quát,một cách t−ơng tự ta có : ∑∑∑ ω=ω== xisinyn2bxicosyn2an y a ii0 Ch−ơng trình tìm các hệ số ai và bi đ−ợc thể hiện nh− sau : Ch−ơng trình 11-7 //xap_xi_sin_cos; #include #include #include #include #define max 11 #define pi 3.15159 void main() { int i,j,m,n,t; float x[max],y[max],a[max],b[max]; char ok; float omg,t1; clrscr(); printf("PHUONG PHAP BINH PHUONG TOI THIEU"); printf("\n"); printf("Cho so so hang sin-cos m = "); scanf("%d",&m); printf("Cho chu ki T = "); scanf("%f",&t1); printf("So diem da cho n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } x[0]=1.0; printf("%4cBANG SO LIEU\n",' '); printf("%8cx%30cy\n",' ',' '); for (i=1;i<=n;i++) 200 printf("%4c%8.4f%23c%8.4f\n",' ',x[i],' ',y[i]); ok=' '; t=1; flushall(); while (t) { printf("Co sua so lieu khong(c/k): "); scanf("%c",&ok); if (toupper(ok)=='C') { printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); flushall(); } if (toupper(ok)!='C') t=0; } printf("\nCAC GIA TRI DA CHO\n"); printf("\n"); printf(" X Y\n"); for (i=1;i<=n;i++) printf("%c%8.3f%c%8.3f\n",' ',x[i],' ',y[i]); printf("\n"); a[0]=0.0; omg=2*pi/t1; for (i=1;i<=n;i++) a[0]+=y[i]; a[0]/=n; for (j=1;j<=m;j++) { a[j]=0.0; for (i=1;i<=n;i++) a[j]+=y[i]*cos(j*omg*x[i]); a[j]=2*a[j]/n; } for (j=1;j<=m;j++) { b[j]=0.0; for (i=1;i<=n;i++) b[j]+=y[i]*sin(j*omg*x[i]); b[j]=2*b[j]/n; } printf("\n"); printf("TAN SO GOC OMEGA = %10.5f\n",omg); 201 printf("HE SO HANG\n"); printf("a[0] = %8.4f\n",a[0]); printf("CAC HE SO BAC CAO\n"); printf("%5ccos%25csin\n",' ',' '); for (i=1;i<=m;i++) printf("%8.4f%21c%6.4f\n",a[i],' ',b[i]); getch(); } Với hàm cho bằng bảng số : x 0 0.15 0.3 0.45 0.6 0.75 0.9 1.05 1.2 1.3 y 2.200 1.595 1.03 1 0.72 2 0.78 6 1.20 0 1.80 5 2.36 9 2.67 8 2.614 Chọn số hệ số sin-cos m = 1,số điểm cho tr−ớc n = 10,chu kì T = 15 ta nhận đ−ợc kết quả tính a0 = 1.7 ; a1 = 0.5 ; b1 = -0.8661 và ω = 4.18879.Nh− vậy hàm xấp xỉ có dạng : f(x) = 1.7 + 0.5cos(4.18879x) - 0.8661sin(4.18879x) 5.Hàm hữu tỉ : Khi quan hệ y = f(x) có dạng đ−ờng cong bão hoà hay dạng arctan,tan v.v ta dùng hàm xấp xỉ là hàm hữu tỉ dạng đơn giản : xb ax y += Lấy nghịch đảo của nó ta có : a 1 x 1 a b y 1 += Đặt 1/y = Y,1/x = X,b/a = B và 1/a = A ph−ơng trình trên sẽ có dạng : Y = A + BX và là một đa thức bậc một.Do vậy ta có hệ ph−ơng trình đối với các hệ số A và B là : ⎪⎪⎩ ⎪⎪⎨ ⎧ =+ =+ ∑ ∑∑ ∑ ∑ = == = = n 1i n 1i ii 2 i n 1i i n 1i n 1i ii yx 1 x 1 B x 1 A y 1 x 1 BnA và từ đó tính đ−ợc a và b.Ch−ơng trình sau mô tả thuật toán trên Ch−ơng trình 11-.8 //xap xi huu_ty; #include #include #include #include #define k 11 void main() { float x[k],y[k]; float a,b,a1,b1,c,d,e; int i,n,t; 202 char ok; clrscr(); printf("PHUONG PHAP BINH PHUONG TOI THIEU"); printf("\n"); printf("So diem da cho n = "); scanf("%d",&n); for (i=1;i<=n;i++) { printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); } printf("%4cBANG SO LIEU\n",' '); printf("%8cx%30cy\n",' ',' '); for (i=1;i<=n;i++) printf("%4c%8.4f%23c%8.4f\n",' ',x[i],' ',y[i]); ok=' '; t=1; flushall(); while (t) { printf("Co sua so lieu khong(c/k): "); scanf("%c",&ok); if (toupper(ok)=='C') { printf("Chi so cua phan tu can sua i = "); scanf("%d",&i); printf("Gia tri moi : "); printf("x[%d] = ",i); scanf("%f",&x[i]); printf("y[%d] = ",i); scanf("%f",&y[i]); flushall(); } if (toupper(ok)!='C') t=0; } printf("CAC GIA TRI DA CHO\n"); printf("\n"); printf("X = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ',x[i]); printf("\n"); printf("Y = "); for (i=1;i<=n;i++) printf("%c%8.3f",' ' ,y[i]); printf("\n"); a=n; 203 b=0.0; c=0.0; d=0.0; e=0.0; for (i=1;i<=n;i++) { b+=1/x[i]; c+=1/y[i]; d+=1/(x[i]*x[i]); e+=1/(x[i]*y[i]); } a1=(c*d-b*e)/(a*d-b*b); b1=(a*e-b*c)/(a*d-b*b); a=1/a1; b=b1*a; printf("\n"); printf("Cac he so cua ham huu ty\n"); printf("a = %10.5f b = %10.5f",a,b); getch(); } Với dãy số liệu đã cho : x 1 2 3 4 5 y 0.333333 3 0.5 0.6 0.66666 0.7142857 ta nhận đ−ợc từ ch−ơng trình trị số a = 1 và b = 2