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

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.