Đề 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
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.
Bạn đang xem nội dung tài liệu Đề 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, để 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
y vẽ = H max- (f(x) - ymin)*xy
if'i = 1 -> move to (i + lmin, y vẽ)
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]) - y min)
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.
