Đề tài Khớp đường cong bằng phương pháp bình phương nhỏ nhất

2.4. Xây dựng hàm tính giá trị

2.4.1. Hàm tính luỹ thừa

Fuction luỹ thừa (x: Real; n: integer): real.

S : Real

i : integer

S = 1

For i = 1 n

S = S*x

Return S.;

2.4.2. Xây dựng hàm tính f(x)

For i = m 0

S = S + C (i + 1) + luỹ thừa (x, i)

 

doc15 trang | Chia sẻ: netpro | Ngày: 06/04/2013 | Lượt xem: 1138 | Lượt tải: 1download
Bạn đang xem nội dung tài liệu Đề tài Khớp đường cong bằng phương pháp bình phương nhỏ nhất, để tải tài liệu về máy bạn click vào nút DOWNLOAD ở trên
Lêi nãi ®Çu Ngµy nay thÕ giíi chóng ta ®· vµ ®ang lµ kû nguyªn cña sù bïng næ th«ng tin. Trong thêi ®¹i ngµy nay khoa häc kü thuËt ngµy mét ®­îc øng dông réng r·i kh¾p mäi n¬i trong mäi lÜnh vùc, mäi nghµnh nghÒ vµ tin häc ph¸t triÓn kh«ng ngõng tõng gi©y, tõng phót.M¸y tÝnh ®iÖn tö ®­¬c sö dông ,khai th¸c triÖt ®Ó trong mäi c«ng viÖc , nh­ trong hµng kh«ng , qu©n sù , ng©n hµng ,®o l­êng , tÝnh to¸n , vµ nã còng ®­îc sö dông ®Ó gi¶i quyÕt c¸c bµi to¸n tõ ®¬n gi¶n ®Õn phøc t¹p.§Ò tµi thùc tËp tèt nghiÖp nµy còng ®ù¬c ¸p dông ®Ó thùc hiÖn , tªn ®Ò tµi: Cµi ®Æt bµi to¸n:"Khíp ®­êng cong" b»ng ph­¬ng ph¸p b×nh ph­¬ng nhá nhÊt. §Ó hoµn thµnh ®Ò tµi tèt nghiÖp nµy, bªn c¹nh nç lùc cña b¶n th©n em cßn nhËn ®­îc sù gióp quý b¸u cña c¸c thÇy, c« vµ c¸c b¹n. Em xin ch©n thµnh c¶m ¬n thÇy D­¬ng Th¨ng Long ®· tËn t×nh h­íng dÉn em trong suèt qu¸ tr×nh thùc hiÖn ®Ò tµi. C¶m ¬n khoa C«ng nghÖ tin häc ®· t¹o mäi ®iÒu kiÖn gióp ®ì em trong suèt qu¸ tr×nh häc tËp ë tr­êng vµ c¸c thÇy c« ®· hÕt lßng truyÒn ®¹t cho em nh÷ng kiÕn thøc vµ kinh nghiÖm quý gi¸. Còng nh­ tÊt c¶ c¸c b¹n ®· cïng trao ®æi gióp ®ì nhau trong qu¸ tr×nh häc tËp. C¬ së vµ c¸c thuËt to¸n trong ®Ò tµi "Khíp ®­êng cong" b»ng ph­¬ng ph¸p b×nh ph­¬ng nhá nhÊt 1. C¬ së: Khi nh÷ng gi¸ trÞ d÷ liÖu kh«ng chÝnh x¸c, th­êng ta cÇn ph¶i t­ëng t­îng ra h×nh d¹ng cña hµm khíp vµ d÷ liÖu. Hµm nµy cã thÓ phô thuéc mét tham sè. f(x) = f(c1,c2...., cM, x) vµ thñ tôc ®iÒu chØnh ®­êng cong lµ t×m c¸ch chän c¸c tham sè nµo khíp nhÊt víi c¸c gi¸ trÞ quan s¸t ®­îc ë c¸c ®iÓm ®· cho. NÕu hµm lµ mét ®a thøc (víi tham sè lµ nh÷ng hÖ sè) vµ c¸c gi¸ trÞ lµ chÝnh x¸c th× ®©y chÝnh lµ phÐp néi suy. Nh­ng ta ®ang xÐt nh÷ng hµm tæng qu¸t h¬n vµ d÷ liÖu kh«ng ®óng. §Ó gi¶n ®¬n ho¸ vÊn ®Ò, ta tËp trung vµo viÖc nèi nh÷ng hµm lµ tæ hîp tuyÕn tÝnh cña hµm ®¬n gi¶n h¬n, víi c¸c tham sè Èn lµ c¸c hÖ sè: f(x) = c1f1(x) + c2f2(x) + ........ +cMfM(x) Hµm nµy bao hµm hÇu hÕt c¸c hµm mµ ta quan t©m. Sau khi nghiªn cøu xong tr­êng hîp nµy, ta sÏ xÐt ®Õn nh÷ng hµm tæng qu¸t h¬n. Mét c¸c phæ th«ng ®Ó ®o møc ®é tèt cña hµm nèi lµ tiªu chuÈn b×nh ph­¬ng tæng qu¸t. ë ®©y sai sè ®­îc tÝnh b»ng c¸ch thªm vµo b×nh ph­¬ng sai sè ë mçi ®iÓm quan s¸t ®­îc. 1£ j £N E = Z (f(xj) - þ)2 §©y lµ mét phÐp ®o rÊt tù nhiªn: b×nh ph­¬ng ®­îc thùc hiÖn ®Ó lµm khö c¸c ­íc l­îc gi÷a c¸c sai sè víi c¸c sè hiÖu kh¸c nhau. HiÓn nhiªn, ®iÒu ng­êi ta mong muèn nhÊt lµ t×m ra mét phÐp chän c¸c tham sè sao cho tèi thiÓu ho¸ E. §iÒu nµy khiÕn phÐp chän cã thÓ ®­îc tÝnh cã hiÖu qña; ®©y lµ ph­¬ng ph¸p b×nh ph­¬ng nhá nhÊt, ph­¬ng ph¸p nµy ®óng nh­ ®Þnh nghÜa cña nã. §Ó ®¬n gi¶n ho¸ ®¹o hµm, xÐt tr­êng hîp M=2, N=3. Gi¶ sö cã ba ®iÓm x1, x2,x3 vµ c¸c gi¸ trÞ t­¬ng øng y1,y2,y3 tho¶ hµm cã d¹ng: f(x) = (c1f1(x1) + c2f2(x1) - y1)2 + (c1f1(x2) + c2f2(x2) - y2)2 + (c1f1(x3) + c2f2(x3) - y3)2 §Ó t×m c¸c phÐp chän cña c1 vµ c2 sao chá tèi tiÓu ho¸ sao sè nµy, ®¬n gian chØ cÇn g¸n zero c¸c ®¹o hµm dE/dc1 vµ dE/dc2 Víi c1 ta cã : dE/dc1 = 2(c1f1(x1) + c2f2(x1) - y1)2 f1(x1) + (c1f1(x2) + c2f2(x2) - y2)2 f2(x2) + (c1f1(x3) + c2f2(x3) - y3)2 f3(x3) ViÖc g¸n ®¹o hµm b»ng zero cho ra ph­¬ng tr×nh cã biÕn lµ c1 vµ c2 ph¶i tho¶ (f1(x1), lµ c¸c "h»ng sè" cã trÞ biÕt tr­íc). c1(f1(x1)f1(x1) + f1(x2)f1(x2) + (f1(x3)f1(x3) + c2(f2(x1)f1(x1) + (f2(x2)f1(x2) + (f2(x3)f1(x3)) = y1f1 (x1) + y2f1 (x2) + y3f1 (x3) Ta cã ph­¬ng tr×nh t­¬ng tù, khi g¸n ®¹o hµm dE/dc2 vÒ zero Gi¶ sö ®a thøc xÊp xØ bËc m cña f lµ P(x) = Sj=0 Cjx víi Cj (j = O1m) x¸c ®Þnh tõ hÖ D0C0 +D1C1 + ...... + DmCm = t0 D1C0 +D2C1 + ...... + Dm+1Cm = t1 .............. DmC0 +Dm+1C1 + ...... + D2mCm = tm Trong ®ã Dk = S xg (k= 0,2m) ti = S xg f (xg) (i = o1m) víi Cj (j = 0,m) x¸c ®Þnh tõ hÖ ph­¬ng tr×nh trªn trong ®ã. 2. C¸c thuËt to¸n dïng trong ch­¬ng tr×nh Cµi ®Æt bµi to¸n: "Khíp ®­êng cong" b»ng ph­¬ng ph¸p b×nh ph­¬ng nhá nhÊt 2.1. NhËp d÷ liÖu + nhËp tõ file: NhËp c¸c to¹ ®é xi , yi vµo m¸y to¹ ®é, sè ®iÓm n, bËc m cña ®a thøc. - hái xem cã ghi vµo file kh«ng nÕu hái nhËp tªn tªn file ghi d÷ liÖu vµo file ®ã. + nhËp tõ file: yªu cÇu nhËp vµo tªn file, ®äc file Read (f,n,m) n¹p c¸c qt tõ file vµo hai m¶ng to¹ ®é do x, to¹ ®é y For i = 1-> n Read (f, to¹ ®é x[i], to¹ ®é y[i]) 2.2. TÝnh m¶ng -TÝnh to¸n c¸c hÖ sè cña m¶ng a[i] [j], b[i] theo hÖ pt (1) 2.3.Dïng ph­¬ng ph¸p gauss vµ thuËt to¸n sytru gi¶i hÖ pt (1) Gauss Input n, (¹i)1n; (bi)1n IER = 0 akk ¹ 0 Pi=aik/akk k= 1, 2,...,n i= k + 1, ..., n j= k + 1, ..., n aþj = piakj bi = bi - pibk Call Sytru IER = 1 Print (xi)1n End ThuËt to¸n gauss (§­a ma trËn vÒ d¹ng tam gi¸c ) no yes thuËt to¸n Sytru (gi¶i hÖ ph­¬ng tr×nh d¹ng tam gi¸c) IRE = 1 Procedure Input n, {bi}1n, aij (1£ i £ j £ n) IER = 0 ann ¹ 0 xn = bn / ann i = n-1, n-2,...,1 S = bi j = i + 1, ..., n S = S - aijxij aij ¹ 0 xi = S/aij Return Procedure Input n, {bi}1n, aij (1£ i £ j £ n) IER = 0 ann ¹ 0 xn = bn / ann i = n-1, n-2,...,1 S = bi j = i + 1, ..., n S = S - aijxij aij ¹ 0 xi = S/aij Return no yes no yes 2.4. X©y dùng hµm tÝnh gi¸ trÞ 2.4.1. Hµm tÝnh luü thõa Fuction luü thõa (x: Real; n: integer): real. S : Real i : integer S = 1 For i = 1 ® n S = S*x Return S.; 2.4.2. X©y dùng hµm tÝnh f(x) For i = m ® 0 S = S + C (i + 1) + luü thõa (x, i) 3. VÏ ®å thÞ 3.1. T×m xmin , xmax , ymin , ymax ,trong toµn ®å thÞ vµ c¸c ®iÓm nhËp vµo x¸c ®Þnh kx , xy Kx = (Lmax - Lmin)/ (xmax - xmin) (Hmax - H min) (ymax - ymin) Ky = vÏ ®å thÞ tõ ®iÓm x1-> x2 , x = x1, i = 0 while x<x2 i = i + 1 Begin x = x1 + x2 i/ kx yvÏ = H max- (f(x) - ymin)*xy if'i = 1 -> move to (i + lmin, yvÏ) ESC line to (i + lmin, yvÏ) + vÏ trôc, to¹ ®é + vÏ c¸c ®iÓm nhËp vµo For i = 1 -> n Begin x = kx + (to¹ ®é x[i]) - ymin) y vÏ = Hmax - (Toado y[i]) - ymin)*ky vÏ ®iÓm (x+lmin, yvÏ, to¹ ®é x[i], to¹ ®é y[i]); 4. Ch­¬ng tr×nh b»ng ph­¬ng ph¸p b×nh ph­¬ng M· nguån uses crt,graph; const Hmin=150; Lmin=200; Hmax=480-Hmin; Lmax=640-Lmin; Type mang2=array[1..20,1..20] of real; mang1=array[1..20] of real; mang_diem=array[1..20] of real; var a:mang2; b:mang1; n,m:integer; toadox,toadoy:mang_diem; c:mang1; ok:boolean; r:real; (*----------------------------*) procedure hien_mang(so:integer;a:mang2;b:mang1); var i,j:integer; Begin {for i:=1 to so do begin for j:=1 to so do write(a[i,j]:0:2,' '); write(' | ',b[i]:0:2); writeln; end; } End; (*----------------------------*) function luythua(x:real;n:integer):real; var i:integer; s:real; Begin s:=1; if n>0 then for i:=1 to n do s:=s*x; luythua:=s; End; (*----------------------------*) procedure nhap_tu_ban_fim(var n,m:integer;var toadox,toadoy:mang_diem); var i:integer; x,y:real; ch:char; tenfile:string; f:text; Begin textcolor(14); textbackground(1); window(1,1,80,24); clrscr; write('Nhap vao so diem:');readln(n); writeln('Nhap vao cac diem:'); for i:=1 to n do begin write('vao toa do cua diem ',i,' x=');readln(toadox[i]); write('vao toa do cua diem ',i,' y=');readln(toadoy[i]); end; write('Nhap vao so bac:');readln(m); Writeln('Nhap Du Lieu OK !'); repeat Write('Ban CO Muon Ghi Du Lieu Vao File Khong (C/K)?:'); ch:=readkey; until upcase(ch) in['C','K']; writeln; if upcase(ch)='C' then Begin write('Moi nhap vao ten file:'); readln(tenfile); assign(f,tenfile); {$I-} rewrite(f); if IOResult 0 then begin gotoxy(30,10); Writeln('File not found'); readln; exit; end; {$I+} writeln(f,n,' ',m); for i:=1 to n do write(f,toadox[i],' ',toadoy[i],' '); close(f); Writeln('Nhap File OK !'); readkey; end; end; (*----------------------------------*) procedure nhap_tu_file(var n,m:integer; var toadox,toadoy:mang_diem); var f:text; tenfile:string; i:integer; Begin n:=0; m:=0; write('Moi nhap vao ten file:'); readln(tenfile); assign(f,tenfile); {$I-} reset(f); if IOResult 0 then begin gotoxy(30,10); Writeln('File not found'); readln; exit; end; {$I+} readln(f,n,m); for i:=1 to n do read(f,toadox[i],toadoy[i]); close(f); Writeln(' Nhap Tu File Ok !'); readkey; End; (*----------------------------------*) procedure tinh_mang; var i,j,dem:integer; so:array[0..200] of real; Begin dem:=0; for i:=0 to 2*m do begin so[i]:=0; for j:=1 to n do so[i]:=so[i]+luythua(toadox[j],i); end; for i:=0 to m do begin dem:=i*2-1; for j:=i to m do begin dem:=dem+1; a[i+1,j+1]:=so[dem]; a[j+1,i+1]:=so[dem]; end; end; for i:=0 to m do begin b[i+1]:=0; for j:=1 to n do b[i+1]:=b[i+1]+toadoy[j]*luythua(toadox[j],i); end; End; (*--------------- Giai Phuong Trinh --------------*) procedure gauss(m:integer;a:mang2;b:mang1; var nghiem:mang1;var ok:boolean); var k,i,j:integer; s:real; mi:mang1; Begin ok:=true; for i:=1 to m+1 do Begin if a[i,i]=0 then begin ok:=false;exit;end; for j:=i+1 to m+1 do mi[j]:=a[j,i]/a[i,i]; for j:=i+1 to m+1 do begin for k:=1 to m+1 do a[j,k]:=a[j,k]-a[i,k]*mi[j]; b[j]:=b[j]-b[i]*mi[j]; end; End; (*----------------Sytru-------------*) nghiem[m+1]:=b[m+1]/a[m+1,m+1]; for i:=m downto 1 do Begin s:=b[i]; for j:=i+1 to m+1 do s:=s-a[i,j]*nghiem[j]; nghiem[i]:=s/a[i,i]; end; End; (*----------------------------------*) function f(x:real;n,m:integer;c:mang1):real; var s:real; i:integer; Begin s:=0; for i:=m downto 0 do s:=s+c[i+1]*luythua(x,i); f:=s; End; (*----------------------------------*) Procedure Tinh_Min_Max_X_Y_(n,m:integer;c:mang1; var kx:real;var xmin,xmax,ymin,ymax:real); var i,j:integer; x:real; Begin xmin:=toadox[1]; xmax:=toadox[1]; ymin:=toadoy[1]; ymax:=toadoy[1]; for i:=1 to n do begin if xmin>toadox[i] then xmin:=toadox[i]; if xmax<toadox[i] then xmax:=toadox[i]; if ymin>toadoy[i] then ymin:=toadoy[i]; if ymax<toadoy[i] then ymax:=toadoy[i]; end; i:=0; if xmaxxmin then begin kx:=abs((lmax-lmin)/(xmax-xmin)); while x<xmax do begin i:=i+1; x:=xmin+(i/kx); if ymin>f(x,n,m,c) then ymin:=f(x,n,m,c); if ymax<f(x,n,m,c) then ymax:=f(x,n,m,c); end; end; End; (*----------------------------------*) procedure draw_point(x,y:integer;x1,y1:real); var i:integer; s1,s2:string; Begin setcolor(13); line(x-2,y-2,x+2,y+2); line(x+2,y-2,x-2,y+2); if(x>320) then i:=10 else i:=-100; Str(x1:0:1,s1); Str(y1:0:1,s2); s1:='('+s1+','+s2+')'; outtextxy(x+i,y,s1); End; (*----------------------------------*) Procedure Ve_do_thi(n,m:integer;c:mang1); Var kx,ky:real; x,yve,xmin,xmax,ymin,ymax:real; i,j:integer; gd,gm:integer; s1,s2:string; Begin Tinh_Min_Max_X_Y_(n,m,c,kx,xmin,xmax,ymin,ymax); gd:=detect; initgraph(gd,gm,''); setfillstyle(1,0); bar(0,0,getmaxx,getmaxy); if yminymax then begin i:=0; ky:=abs((hmax-hmin)/(ymax-ymin)); setcolor(13); x:=xmin; if xminxmax then begin (* Ve do Thi *) x:=xmin; setcolor(14); i:=0; while x<xmax do Begin x:=xmin+i/kx; yve:=hmax-(f(x,n,m,c)-ymin)*ky; if i=0 then moveto(lmin,round(yve)); lineto(lmin+i,round(yve)); inc(i); end; for i:=1 to n do Begin x:=kx*(toadox[i]-xmin); yve:=hmax-(toadoy[i]-ymin)*ky; draw_point(round(x+lmin),round(yve),toadox[i],toadoy[i]); end; if round((-xmin)*kx)+lmin<640 then if round((-xmin)*kx)+lmin>=0 then begin setcolor(15); yve:=hmax-(ymax-ymin)*ky; line(round((-xmin)*kx)+lmin,5,round((-xmin)*kx+lmin),480); line(round((-xmin)*kx)+lmin,5,round((-xmin)*kx)+lmin-3,5+3); line(round((-xmin)*kx)+lmin,5,round((-xmin)*kx)+lmin+3,5+3); setcolor(13); outtextxy(round((-xmin)*kx)+lmin+10,20,'Y'); end; yve:=hmax+(ymin)*ky; (* Ve truc Hoanh *) if (yve>0)and(yve<480) then begin setcolor(15); line(5,round(yve),639,round(yve)); line(639,round(yve),639-3,round(yve)+3); line(639,round(yve),639-3,round(yve)-3); setcolor(13); outtextxy(600,round(yve)-20,'X'); end; End else line(320,hmin,320,hmax); end else Begin if xminxmax then line(lmin+10,240,lmax-10,240) else putpixel(320,240,14); End; readkey; closegraph; End; (*----------------------------------*) procedure Menu; var ch:char; i:integer; Begin gotoxy(10,3); Writeln('Chuong Trinh Khop Ham so Bang Phuong Phap Binh Phuong Nho Nhat'); gotoxy(10,4); Writeln('---------------------------- =o0o= ----------------------------'); gotoxy(30,5); Writeln(' 1. Nhap Du Lieu Tu Ban Phim '); gotoxy(30,6); Writeln(' 2. Nhap Tu File '); gotoxy(30,7); Writeln(' 3. Xem Cac He So '); gotoxy(30,8); Writeln(' 4. Ve Do Thi '); gotoxy(30,9); Writeln(' 5. Thoat Khoi Chuong Trinh '); gotoxy(25,11); Write('Ban Hay Chon Mot Trong Cac So:'); repeat ch:=readkey; case ch of '1':begin clrscr; nhap_tu_ban_fim(n,m,toadox,toadoy); Writeln('Nap Du Lieu Ok !'); tinh_mang; gauss(m,a,b,c,ok); Readln; exit; end; '2':begin clrscr; nhap_tu_file(n,m,toadox,toadoy); Writeln('Nap Du Lieu Ok !'); tinh_mang; gauss(m,a,b,c,ok); Readln; exit; end; '3':begin clrscr; if n>0 then begin hien_mang(m+1,a,b); for i:=m+1 downto 1 do writeln('c[',i,']:',c[i]:0:4); end else Writeln('Ban Chua Vao Du Lieu !'); readkey; exit; end; '4':begin clrscr; if (n>0)then begin if ok then ve_do_thi(n,m,c) else Writeln('Khong the khop ham'); end else Writeln('Ban Chua Nap Du Lieu'); readkey; exit; end; '5':halt(1); end; until false; End; (*------------ Chuong Trinh Chinh------------*) Begin n:=0; m:=0; repeat textcolor(14); textbackground(1); clrscr; menu; until false; End.

Các file đính kèm theo tài liệu này:

  • docCài đặt bài toán- Khớp đường cong bằng phương pháp bình phương nhỏ nhất.DOC