Giáo trình phương pháp số

t;=m;i++) for(j=0;j<i;j++) aa[i][j]=aa[j][i]; //Tinh ve phai cua he phuong trinh for(i=0;i<=m;i++) {aa[i][m+1]=ft[i][0]*yqs[0]; for(k=1;k<=nqs;k++) aa[i][m+1]+=ft[i][k]*yqs[k]; } for(i=0;i<=nqs;i++) delete [] f[i]; for(i=0;i<=m;i++) delete [] ft[i]; gjordan(aa,a,m); //Tinh tong binh phuong cac sai so double ss,fa,xx; ss=0; for(i=0;i<=nqs;i++) {fa=1;xx=1; for(j=1;j<=m;j++) {xx=xx*xqs[i];fa+=a[j]*xx;} ss+=(yqs[i]-fa)*(yqs[i]-fa); } return ss; } //=============================================== void hoiquy() {int i,j,mtmp;double h;char ordst[2],mchon; if(nqs<2) {cout<<"\nCo le ban chua nhap so lieu!";delay(2000);return;} while(true) {clrscr(); cout>mtmp; if(mtmp<1) {cout<<"Bac da thuc phai khong am";delay(1000);continue;} if(mtmp>17) {cout<<"Bac da thuc qua lon";delay(1000);continue;} clrscr(); break; } mord=mtmp; //Ve do thi double scalex,scaley; /*[a,b] la khoang chua cac gia tri x, scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx,gocy,n,color; /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dijnh tuy thuoc vao dang ham*/ double (*f)(double);//Ham can ve do thi double a,b,ss;char ssst[]=""; int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); while(true) {ss=regresspoli(apoli,mord); sprintf(ssst,"%8.2f",ss); for(i=mord;i>=0;i--) if(apoli[i]!=0) break; //Xac dinh bac cua da thuc hoi quy itoa(i,ordst,10);//Chuyen sang dang ky tu char ftitle[]="Da thuc hoi quy bac "; strcat(ftitle,ordst);//ftitle se la dong chu co them bac cua da thuc cleardevice(); setbkcolor(GREEN); setcolor(BLUE); gocx=gocy=0;//De chhuong trinh tu xac dinh n=100;//100 diem chia f=poli; a=b=xqs[0];//Tim min va max cua xqs(i) for(i=1;ib?xqs[i]:b;} //Ve cac diem quan sat vequansatxy(xqs,yqs,nqs,gocx,gocy,scalex,scaley,RED,2); //Ve do thi da thuc hoi quy vedothi(f,a,b,n,gocx,gocy,scalex,scaley,RED,0); setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,"PHUONG PHAP BINH PHUONG CUC TIEU"); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve chu thich ham va da thuc noi suy settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); setcolor(RED); int x1,y1,x2,y2; x1=(getmaxx()-gocx)/2-25; y1=30-gocy; x2=x1+50; y2=y1; line(x1,y1,x2,y2); line(x1,y1+1,x2,y2+1);//Ve cho dam hon outtextxy(x2+10,y2,"Da thuc hoi quy"); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Sai so:"); outtextxy(x2+10,y2,ssst); setcolor(BLUE); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"PgUp:"); outtextxy(x2+10,y2,"Tang 1 bac"); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"PgDn:"); outtextxy(x2+10,y2,"Giam 1 bac"); //in anh //inanh(1,gocx-getmaxx(),gocy-getmaxy(),getmaxx()-gocx,getmaxy()-gocy); mchon=getch(); if(mchon==0) mchon=getch(); if(mchon==73 && mord<17)//Bam PgUp {mord++;continue;} if(mchon==81 && mord>0) //Bam PgDn {mord--;continue;} break; } closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"XAC DINH DA THUC HOI QUY BANG PHUONG PHAP BPCT"); fprintf(fp,"\nSo mau quan sat n = %d",nqs); fprintf(fp,"\n i x[i] y[i]"); for(i=0;i<=nqs;i++) fprintf(fp,"\n%2d %8.2f %8.2f",i,xqs[i],yqs[i]); fprintf(fp,"\n\nKET QUA TINH TOAN CAC PHUONG PHAP BPCT:"); fprintf(fp,"\n========================================"); fprintf(fp,"\nBac cua da thuc hoi quy m = %d",mord); fprintf(fp,"\nCac he so a[0],a[1],...,a[m]: "); for(i=0;i<=mord;i++) fprintf(fp,"%8.2f ",apoli[i]); fprintf(fp,"\nGIA TRI CAC DA THUC HOI QUI TAI MOT SO DIEM:"); fprintf(fp,"\n x[i] y[i] p(x[i])"); for(i=0;i<=nqs;i++) fprintf(fp,"\n%8.2f %8.2f %8.2f",xqs[i],yqs[i],poli(xqs[i])); fclose(fp); xemtep(resultfile); gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; } //=============================================== void hoiquyhamdb() {double a,b,h;int i,j,mm; while(true) {textbackground(BLUE); textcolor(YELLOW); clrscr(); gotoxy(2,2); cout<<endl<<" 1. Ham sin"; cout<<endl<<" 2. Ham sin(x)*sin(x)*sin(x)"; cout<<endl<<" 3. Ham mat do phan bo chuan"; cout<<endl<<" z. Ket thuc"; cout z de chon"; char mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr();//So lieu co dang ham sin a=-2*M_PI;b=2*M_PI; nqs=40;mm=20; h=(b-a)/nqs; for(i=0;i<=nqs;i++) {xqs[i]=a+i*h;yqs[i]=sin(xqs[i%20]);} randomize(); for(i=0;i<=nqs;i++) yqs[i]+=double(random(20))/80-0.125; gotoxy(38,25);cout<<"Da tao xong, nhan phim bat ky de tiep tuc"; getch(); return; case '2':clrscr();//So lieu co dang ham sin mu 3 a=-2*M_PI;b=2*M_PI; nqs=40;mm=20; h=(b-a)/nqs; for(i=0;i<=nqs;i++) {xqs[i]=a+i*h;yqs[i]=sinmu3(xqs[i%20]);} randomize(); for(i=0;i<=nqs;i++) yqs[i]+=double(random(20))/80-0.125; gotoxy(38,25);cout<<"Da tao xong, nhan phim bat ky de tiep tuc"; getch(); return; case '3':clrscr();//Noi suy ham mat do phan bo chuan a=gausstb-3*gaussdlc;b=gausstb+3*gaussdlc; nqs=40;mm=20; h=(b-a)/nqs; for(i=0;i<=nqs;i++) {xqs[i]=a+i*h;yqs[i]=gauss(xqs[i%20]);} randomize(); for(i=0;i<=nqs;i++) yqs[i]+=double(random(20))/80-0.125; gotoxy(38,25);cout<<"Da tao xong, nhan phim bat ky de tiep tuc"; getch(); return; } textbackground(BLACK); textcolor(WHITE); clrscr(); } } //=============================================== void main() {int i,j;char mchon; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(20,2);cout<<"DA THUC HOI QUY"; gotoxy(2,2); cout<<endl<<" 1. Nhap so lieu x,y"; cout<<endl<<" 2. Xem so lieu"; cout<<endl<<" 3. Thuc hien hoi quy "; cout<<endl<<" 4. Tao so lieu tu mot so ham quen thuoc"; cout<<endl<<" z. Ket thuc"; cout z de chon"; fflush(stdin); mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr(); nhapsl(xqs,yqs,nqs); break; case '2':clrscr(); xemsl(xqs,yqs,nqs); break; case '3':clrscr(); hoiquy(); break; case '4':clrscr(); hoiquyhamdb(); break; } //gotoxy(50,25); //cout<<"Nhan phim bat ky de tiep tuc"; //getch(); } } //=============================================== /*Giai he phuong trinh tuyen tinh bang cach dua ma tran a ve dang ma tran don vi. Tra ve gia tri true neu co nghiem */ int gjordan(kmatran a,double *x,int n) {int i,j,k,h;double tmp,p; int n1=n+1; for(i=0;i<=n;i++) //Vong lap cac buoc khu {//Tim hang co phan tu dau lon nhat h=i; for(k=i+1;k<=n;k++) if(fabs(a[k][i])> fabs(a[h][i])) h=k; if(a[h][i]==0) {cout<<"Ma tran suy bien";delay(1000);return false;} if(h!=i) //Doi hang i va hang h vi a[h][i] > a[i][i] {int j;double tmp; for(j=i;j<=n1;j++) {tmp=a[i][j];a[i][j]=a[h][j];a[h][j]=tmp;} } //chuyen he so a[i][i] = 1 tmp=a[i][i]; for(j=i;j<=n1;j++) a[i][j] = a[i][j]/tmp; //Bat tinh lai cac hang for(k=0;k<=n;k++) {if(k==i) continue; p=a[k][i]; /*Vi ta biet a[k][i] =0 sau bien doi, chi tinh tu a[k][i+1]*/ for(j=i+1;j<=n1;j++) a[k][j]=a[k][j] - p*a[i][j]; } } for(i=0;i<=n;i++) x[i]=a[i][n+1]; return true; } //=============================================== double sinmu3(double x) {return sin(x)*sin(x)*sin(x);} //=============================================== //a=-0.5;b=2; double gauss(double x) {return (1/(gaussdlc*sqrt(2*M_PI)))*(exp(-(x-gausstb)*(x-gausstb)/ (2*gaussdlc*gaussdlc))); } Tính gần đúng nghiệm của phương trình phi tuyến Phương pháp chia đôi tìm nghiệm: CHIADOI.CPP //CHIADOI.CPP /*Phuong phap tim nghiem bang cach chia doi lien tiep khoang [a,b]*/ #include "zheader.cpp" #include "zlibrary.cpp" char *ftitle,*datafile="CHIADOI.DAT",*resultfile="CHIADOI.KQ"; typedef double kvecto[200]; double epsix=1.0E-03,epsiy=1.0E-03; int kmax=20; //a=-0.5;b=2; double g(double x) {return sin(x)-x*x*cos(x); } //a=1;b=2; double h(double x) {return x*x*x-x-1; } //=============================================== //a=-1.2;b=2; double p(double x) {return x*x-2; } //=============================================== int chiadoi(double (*f)(double),double a,double b,double &x, double &errx,double &erry,int &buoclap) /*Phuong phap chia doi de tim nghiem f(x)=0 trong khoang trai dau [a,b]. Neu co nghiem thi tra ve gia tri true, neu khong thi tra ve gia tri false. x la nghiem, errx, erry sai so tren x va tren y, buoclap la so buoc lap da thuc hien*/ {clrscr(); {int c=2,h=2; gotoxy(c,h); printf("CAC THONG SO DE TINH TOAN:"); gotoxy(c,++h); printf("Can duoi a = : %6.1f",a); gotoxy(c,++h); printf("Can tren b = : %6.1f",b); gotoxy(c,++h); printf("So buoc lap toi da: %2d",kmax); gotoxy(c,++h); printf("Sai so |x[n]-x[0]|: %6.5f",epsix); gotoxy(c,++h); printf("Sai so |f(x[n])| : %6.5f",epsiy); gotoxy(c,h+2); printf("Co chap nhan khong (C/K)"); char mchon=toupper(getch()); if(mchon=='K') {h+=2; gotoxy(c,h); printf("NHAP LAI CAC THONG SO: "); gotoxy(c,++h); printf("Can duoi a = :");scanf("%lf",&a); gotoxy(c,++h); printf("Can tren b = :");scanf("%lf",&b); gotoxy(c,++h); printf("So buoc lap toi da:");scanf("%d",&kmax); gotoxy(c,++h); printf("Sai so |x[n]-x[0]|:");scanf("%lf",&epsix); gotoxy(c,++h); printf("Sai so |f(x[n])| :");scanf("%lf",&epsiy); } } clrscr(); if(f(a)==0) {x=a;errx=erry=0;buoclap=0;return true;} if(f(b)==0) {x=b;errx=erry=0;buoclap=0;return true;} if(f(a)*f(b)>0) {cout<<endl<<"f(a) va f(b) khong trai dau";delay(1000);return false;} int k=1;double c,aa,bb;aa=a;bb=b; kvecto xc; xc[0]=a;xc[1]=b; while(true) {c = (a+b)/2; if(f(a)*f(c)<0) b=c;else a=c; xc[++k]=c; if(b-a<epsix && fabs(f(c))<epsiy) break;//f(c) = 0 if(k>kmax) {cout<<endl<<"Lap chua hoi tu sau "<<kmax<<" buoc"; delay(1000);return false; } } x=c; errx=b-a;erry=fabs(f(c));buoclap=k; //Ve do thi double scalex,scaley; /*scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx=0,gocy=0,n=100,color=RED;//100 diem chia /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dinh tuy thuoc vao dang ham*/ int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); cleardevice(); setbkcolor(GREEN); setcolor(BLUE); vedothi(f,aa,bb,n,gocx,gocy,scalex,scaley,RED,1); //Ve cac diem int xx,yy; color=MAGENTA; setcolor(color); setfillstyle(SOLID_FILL,YELLOW); for(int i=0;i<=buoclap;i++) {xx=xc[i]*scalex; yy=-f(xc[i])*scaley; circle(xx,0,2); floodfill(xx,0,color); circle(xx,yy,2); floodfill(xx,yy,color); line(xx,0,xx,yy); delay(5); } //Ve cac duong tung do mau trang color=WHITE; setcolor(color); xx=a*scalex; yy=-f(a)*scaley; line(xx,0,xx,yy); xx=b*scalex; yy=-f(b)*scaley; line(xx,0,xx,yy); for(i=0;i<=buoclap;i++) {xx=xc[i]*scalex; yy=-f(xc[i])*scaley; line(xx,0,xx,yy); delay(5); } //Ve cac thong bao ket qua setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,"NGHIEM CUA PHUONG TRINH PHI TUYEN"); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve chu thich ham settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); char st1[20],st2[10],st3[10],st4[10],st5[10]; sprintf(st1,"[%4.1f,%4.1f]",aa,bb); sprintf(st2,"%6.4f",x); sprintf(st3,"%2d",buoclap); sprintf(st4,"%6.4f",errx); sprintf(st5,"%6.4f",erry); int x1,y1,x2,y2; x1=(getmaxx()-gocx)/2-25; y1=30-gocy; x2=x1+65; y2=y1; outtextxy(x1,y1,"[a,b] :"); outtextxy(x2+10,y2,st1); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Nghiem x:"); outtextxy(x2+10,y2,st2); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Buoc lap:"); outtextxy(x2+10,y2,st3); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"x-x[0] :"); outtextxy(x2+10,y2,st4); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"f(x) :"); outtextxy(x2+10,y2,st5); getch(); closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"TIM NGHIEM BANG PHUONG PHAP CHIA DOI"); fprintf(fp,"\n(%s)",ftitle); fprintf(fp,"\n\nKET QUA TINH TOAN:"); fprintf(fp,"\n==================="); fprintf(fp,"\nTim nghiem tren khoang [%4.1f,%4.1f]",aa,bb); fprintf(fp,"\nNghiem x[n]: %6.4f",x); fprintf(fp,"\nSo buoc lap: %d",buoclap); fprintf(fp,"\nSai so giua nghiem that va nghiem xap xi: %6.4f",errx); fprintf(fp,"\nGia tri cua f(x[n]): %6.4f",erry); fclose(fp); xemtep(resultfile); return true; } //=============================================== void main() {int buoclap;double a,b,x,errx,erry; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(2,2);cout<<"PHUONG PHAP CHIA DOI GIAI PHUONG TRINH PHI TUYEN"; cout<<endl; cout<<endl<<" 1. sin(x)-x*x*cos(x) = 0";//a=-0.5,b=2; cout<<endl<<" 2. x*x*x-x-1 = 0";//a=1,b=2; cout<<endl<<" 3. x*x-2 = 0";//a=-1,b=2; cout<<endl<<" z. Khong tinh nua"; cout z de chon"; char mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr(); a=-0.5;b=2; ftitle="Tim nghiem pt sin(x)-x*x*cos(x) = 0"; chiadoi(g,a,b,x,errx,erry,buoclap); break; case '2':clrscr(); a=1;b=2; ftitle="Tim nghiem pt x*x*x-x-1 = 0"; chiadoi(h,a,b,x,errx,erry,buoclap); break; case '3':clrscr(); a=-1.2;b=2; ftitle="Tim nghiem pt x*x-2 = 0"; chiadoi(p,a,b,x,errx,erry,buoclap); break; } gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); } } Phương pháp dây cung: DAYCUNG.CPP //DAYCUNG.CPP /*Phuong phap tim nghiem bang phuong phap day cung tren khoang [a,b]*/ #include "zheader.cpp" #include "zlibrary.cpp" char *ftitle,*datafile="DAYCUNG.DAT",*resultfile="DAYCUNG.KQ"; typedef double kvecto[200]; double epsix=1.0E-03,epsiy=1.0E-03; int kmax=20; //a=-0.5;b=2; double g(double x) {return sin(x)-x*x*cos(x); } //a=1;b=2; double h(double x) {return x*x*x-x-1; } //=============================================== //a=-1.2;b=2; double p(double x) {return x*x-2; } //=============================================== int daycung(double (*f)(double),double a,double b,double &x, double &errx,double &erry,int &buoclap) /*Phuong phap day cung de tim nghiem f(x)=0 trong khoang trai dau [a,b]. Neu co nghiem thi tra ve gia tri true, neu khong thi tra ve gia tri false. x la nghiem, errx, erry sai so tren x va tren y, buoclap la so buoc lap da thuc hien*/ {clrscr(); {int c=2,h=2; gotoxy(c,h); printf("CAC THONG SO DE TINH TOAN:"); gotoxy(c,++h); printf("Can duoi a = : %6.1f",a); gotoxy(c,++h); printf("Can tren b = : %6.1f",b); gotoxy(c,++h); printf("So buoc lap toi da: %2d",kmax); gotoxy(c,++h); printf("Sai so |x[n]-x[0]|: %6.5f",epsix); gotoxy(c,++h); printf("Sai so |f(x[n])| : %6.5f",epsiy); gotoxy(c,h+2); printf("Co chap nhan khong (C/K)"); char mchon=toupper(getch()); if(mchon=='K') {h+=2; gotoxy(c,h); printf("NHAP LAI CAC THONG SO: "); gotoxy(c,++h); printf("Can duoi a = :");scanf("%lf",&a); gotoxy(c,++h); printf("Can tren b = :");scanf("%lf",&b); gotoxy(c,++h); printf("So buoc lap toi da:");scanf("%d",&kmax); gotoxy(c,++h); printf("Sai so |x[n]-x[0]|:");scanf("%lf",&epsix); gotoxy(c,++h); printf("Sai so |f(x[n])| :");scanf("%lf",&epsiy); } } clrscr(); if(f(a)==0) {x=a;errx=erry=0;buoclap=0;return true;} if(f(b)==0) {x=b;errx=erry=0;buoclap=0;return true;} if(f(a)*f(b)>0) {cout<<endl<<"f(a) va f(b) khong trai dau";delay(1000);return false;} int k=1;double c,cp,aa,bb;aa=a;bb=b;//cp=c previouse kvecto xa,xb; xa[0]=a;xb[0]=b; c = (a*f(b)-b*f(a))/(f(b)-f(a)); if(f(a)*f(c)<0) b=c;else a=c; xa[1]=a;xb[1]=b; while(true) {cp=c; c = (a*f(b)-b*f(a))/(f(b)-f(a)); if(f(a)*f(c)<0) b=c;else a=c; k++; xa[k]=a;xb[k]=b; if(fabs(c-cp)<epsix && fabs(f(c))<epsiy) break; if(k>kmax) {cout<<endl<<"Lap chua hoi tu sau "<<kmax<<" buoc"; delay(1000);return false; } } x=c; errx=fabs(c-cp);erry=fabs(f(c));buoclap=k; //Ve do thi double scalex,scaley; /*scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx=0,gocy=0,n=100,color=RED;//100 diem chia /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dinh tuy thuoc vao dang ham*/ int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); cleardevice(); setbkcolor(GREEN); setcolor(BLUE); vedothi(f,aa,bb,n,gocx,gocy,scalex,scaley,RED,1); //Ve cac diem va day cung int i,x1,y1,x2,y2; color=MAGENTA; setcolor(color); setfillstyle(SOLID_FILL,YELLOW); for(i=0;i<=buoclap;i++) {x1=xa[i]*scalex; y1=-f(xa[i])*scaley; x2=xb[i]*scalex; y2=-f(xb[i])*scaley; circle(x1,0,2); floodfill(x1,0,color); circle(x1,y1,2); floodfill(x1,y1,color); circle(x2,0,2); floodfill(x2,0,color); circle(x2,y2,2); floodfill(x2,y2,color); line(x1,y1,x2,y2); delay(5); } //Ve cac duong tung do mau trang color=WHITE; setcolor(color); x1=xb[0]*scalex; y1=-f(xb[0])*scaley; line(x1,0,x1,y1); for(i=0;i<=buoclap;i++) {x1=xa[i]*scalex; y1=-f(xa[i])*scaley; line(x1,0,x1,y1); delay(5); } //Ve cac thong bao ket qua setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,"NGHIEM CUA PHUONG TRINH PHI TUYEN"); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve chu thich ham settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); //char *st1,*st2,*st3,*st4,*st5; char st1[20],st2[10],st3[10],st4[10],st5[10]; sprintf(st1,"[%4.1f,%4.1f]",aa,bb); sprintf(st2,"%6.4f",x); sprintf(st3,"%2d",buoclap); sprintf(st4,"%6.4f",errx); sprintf(st5,"%6.4f",erry); //int x1,y1,x2,y2; x1=(getmaxx()-gocx)/2-25; y1=30-gocy; x2=x1+65; y2=y1; outtextxy(x1,y1,"[a,b] :"); outtextxy(x2+10,y2,st1); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Nghiem x:"); outtextxy(x2+10,y2,st2); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Buoc lap:"); outtextxy(x2+10,y2,st3); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"x-x[0] :"); outtextxy(x2+10,y2,st4); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"f(x) :"); outtextxy(x2+10,y2,st5); getch(); closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"TIM NGHIEM BANG PHUONG PHAP DAY CUNG"); fprintf(fp,"\n(%s)",ftitle); fprintf(fp,"\n\nKET QUA TINH TOAN:"); fprintf(fp,"\n==================="); fprintf(fp,"\nTim nghiem tren khoang [%4.1f,%4.1f]",aa,bb); fprintf(fp,"\nNghiem x[n]: %6.4f",x); fprintf(fp,"\nSo buoc lap: %d",buoclap); fprintf(fp,"\nSai so giua nghiem that va nghiem xap xi: %6.4f",errx); fprintf(fp,"\nGia tri cua f(x[n]): %6.4f",erry); fclose(fp); xemtep(resultfile); return true; } //=============================================== void main() {int buoclap;double a,b,x,errx,erry; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(2,2);cout<<"PHUONG PHAP DAY CUNG GIAI PHUONG TRINH PHI TUYEN"; cout<<endl; cout<<endl<<" 1. sin(x)-x*x*cos(x) = 0";//a=-0.5,b=2; cout<<endl<<" 2. x*x*x-x-1 = 0";//a=1,b=2; cout<<endl<<" 3. x*x-2 = 0";//a=-1,b=2; cout<<endl<<" z. Khong tinh nua"; cout z de chon"; char mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr(); a=-0.5;b=2; ftitle="Tim nghiem pt sin(x)-x*x*cos(x) = 0"; daycung(g,a,b,x,errx,erry,buoclap); break; case '2':clrscr(); a=1;b=2; ftitle="Tim nghiem pt x*x*x-x-1 = 0"; daycung(h,a,b,x,errx,erry,buoclap); break; case '3':clrscr(); a=-1.2;b=2; ftitle="Tim nghiem pt x*x-2 = 0"; daycung(p,a,b,x,errx,erry,buoclap); break; } gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); } } Phương pháp Newton: TTUYEN.CPP //TTUYEN.CPP /*Phuong phap tim nghiem bang phuong phap day cung tren khoang [a,b]*/ #include "zheader.cpp" #include "zlibrary.cpp" char *ftitle,*datafile="TTUYEN.DAT",*resultfile="TTUYEN.KQ"; typedef double kvecto[200]; double epsix=1.0E-03,epsiy=1.0E-03; double epsiy1=1.0E-04; //f'(x) phai > epsiy1 int kmax=20; //a=-0.5;b=2; double g(double x) {return sin(x)-x*x*cos(x); } //=============================================== double g1(double x) {return cos(x)-2*x*cos(x)+x*x*sin(x); } //=============================================== //a=1;b=2; double h(double x) {return x*x*x-x-1; } //=============================================== double h1(double x) {return 3*x*x-1; } //=============================================== //a=-1.2;b=2; double p(double x) {return x*x-2; } //=============================================== double p1(double x) {return 2*x; } //=============================================== int ttuyen(double (*f)(double),double (*f1)(double),double a,double b, double &x,double &errx,double &erry,int &buoclap) /*Phuong phap chia doi de tim nghiem f(x)=0 trong khoang trai dau [a,b]. Neu co nghiem thi tra ve gia tri true, neu khong thi tra ve gia tri false. x la nghiem, errx, erry sai so tren x va tren y, buoclap la so buoc lap da thuc hien*/ {clrscr(); {int c=2,h=2; gotoxy(c,h); printf("CAC THONG SO DE TINH TOAN:"); gotoxy(c,++h); printf("Can duoi a = : %6.1f",a); gotoxy(c,++h); printf("Can tren b = : %6.1f",b); gotoxy(c,++h); printf("So buoc lap toi da: %2d",kmax); gotoxy(c,++h); printf("Sai so |x[n]-x[0]|: %6.5f",epsix); gotoxy(c,++h); printf("Sai so |f(x[n])| : %6.5f",epsiy); gotoxy(c,h+2); printf("Co chap nhan khong (C/K)"); char mchon=toupper(getch()); if(mchon=='K') {h+=2; gotoxy(c,h); printf("NHAP LAI CAC THONG SO: "); gotoxy(c,++h); printf("Can duoi a = :");scanf("%lf",&a); gotoxy(c,++h); printf("Can tren b = :");scanf("%lf",&b); gotoxy(c,++h); printf("So buoc lap toi da:");scanf("%d",&kmax); gotoxy(c,++h); printf("Sai so |x[n]-x[0]|:");scanf("%lf",&epsix); gotoxy(c,++h); printf("Sai so |f(x[n])| :");scanf("%lf",&epsiy); } } clrscr(); if(f(a)==0) {x=a;errx=erry=0;buoclap=0;return true;} if(f(b)==0) {x=b;errx=erry=0;buoclap=0;return true;} if(f(a)*f(b)>0) {cout<<endl<<"f(a) va f(b) khong trai dau";delay(1000);return false;} int k=0; kvecto xb;double xp,xn;//previouse x; xb[0]=b; xn=b; while(true) {if(f1(xn)<epsiy1) {cout<<endl<<"Dao ham qua be";delay(1000);return false;} xp=xn; xn = xn-f(xn)/f1(xn); k++; xb[k]=xn; if(fabs(xn-xp)<epsix && fabs(f(xn))<epsiy) break; if(k>kmax) {cout<<endl<<"Lap chua hoi tu sau "<<kmax<<" buoc"; delay(1000);return false; } } errx=fabs(xn-xp);erry=fabs(f(xn));buoclap=k; x=xn; //Ve do thi double scalex,scaley; /*scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx=0,gocy=0,n=100,color=RED;//100 diem chia /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dinh tuy thuoc vao dang ham*/ int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); cleardevice(); setbkcolor(GREEN); setcolor(BLUE); vedothi(f,a,b,n,gocx,gocy,scalex,scaley,RED,1); //Ve cac diem va tiep tuyen int i,x1,y1,x2,y2; color=MAGENTA; setcolor(color); setfillstyle(SOLID_FILL,YELLOW); for(i=0;i<buoclap;i++) {x1=xb[i]*scalex; y1=-f(xb[i])*scaley; x2=xb[i+1]*scalex; y2=-f(xb[i+1])*scaley; circle(x1,0,2); floodfill(x1,0,color); circle(x1,y1,2); floodfill(x1,y1,color); circle(x2,0,2); floodfill(x2,0,color); line(x1,y1,x2,0); delay(5); } //Ve cac duong tung do mau trang color=WHITE; setcolor(color); x1=b*scalex; y1=-f(b)*scaley; line(x1,0,x1,y1); for(i=0;i<=buoclap;i++) {x1=xb[i]*scalex; y1=-f(xb[i])*scaley; line(x1,0,x1,y1); delay(5); } //Ve cac thong bao ket qua setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,"NGHIEM CUA PHUONG TRINH PHI TUYEN"); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve chu thich ham settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); char st1[20],st2[10],st3[10],st4[10],st5[10]; sprintf(st1,"[%4.1f,%4.1f]",a,b); sprintf(st2,"%6.4f",x); sprintf(st3,"%2d",buoclap); sprintf(st4,"%6.4f",errx); sprintf(st5,"%6.4f",erry); x1=(getmaxx()-gocx)/2-25; y1=30-gocy; x2=x1+65; y2=y1; outtextxy(x1,y1,"[a,b] :"); outtextxy(x2+10,y2,st1); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Nghiem x:"); outtextxy(x2+10,y2,st2); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Buoc lap:"); outtextxy(x2+10,y2,st3); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"x-x[0] :"); outtextxy(x2+10,y2,st4); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"f(x) :"); outtextxy(x2+10,y2,st5); getch(); closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"TIM NGHIEM BANG PHUONG PHAP TIEP TUYEN"); fprintf(fp,"\n(%s)",ftitle); fprintf(fp,"\n\nKET QUA TINH TOAN:"); fprintf(fp,"\n==================="); fprintf(fp,"\nTim nghiem tren khoang [%4.1f,%4.1f]",a,b); fprintf(fp,"\nNghiem x[n]: %6.4f",x); fprintf(fp,"\nSo buoc lap: %d",buoclap); fprintf(fp,"\nSai so giua nghiem that va nghiem xap xi: %6.4f",errx); fprintf(fp,"\nGia tri cua f(x[n]): %6.4f\n",erry); fclose(fp); xemtep(resultfile); return true; } //=============================================== void main() {int buoclap;double a,b,x,errx,erry; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(2,2);cout<<"PHUONG PHAP TIEP TUYEN GIAI PHUONG TRINH PHI TUYEN"; cout<<endl; cout<<endl<<" 1. sin(x)-x*x*cos(x) = 0";//a=-0.5,b=2; cout<<endl<<" 2. x*x*x-x-1 = 0";//a=1,b=2; cout<<endl<<" 3. x*x-2 = 0";//a=-1,b=2; cout<<endl<<" z. Khong tinh nua"; cout z de chon"; char mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr(); a=-0.5;b=2; ftitle="Tim nghiem pt sin(x)-x*x*cos(x) = 0"; ttuyen(g,g1,a,b,x,errx,erry,buoclap); break; case '2':clrscr(); a=1;b=2; ftitle="Tim nghiem pt x*x*x-x-1 = 0"; ttuyen(h,h1,a,b,x,errx,erry,buoclap); break; case '3':clrscr(); a=-1.2;b=2; ftitle="Tim nghiem pt x*x-2 = 0"; ttuyen(p,p1,a,b,x,errx,erry,buoclap); break; } gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); } } Tính gần đúng đạo hàm và tích phân xác định Tính xấp xỉ đạo hàm: DAOHAM.CPP //DAOHAM.CPP /*Noi suy dung da thuc Vandermon roi tinh dao ham*/ #include "zheader.cpp" #include "zlibrary.cpp" double gausstb=0, gaussdlc=1; char *datafile="DAOHAM.DAT",*resultfile="DAOHAM.KQ"; typedef double kmatran[30][30],kvecto[50]; kvecto xqs,yqs,avan;int nqs; /*xqs[i] la cac diem quan sat,yqs[i] la ca gia tri quan sat, i=0,1,...,nqs nqs la bac cua da thuc can xac dinh.Khai bao la bien toan cuc de co the thuc hien mot so thao tac chung. anew la cac he so cua da thuc Newton, avan la cac he so tu Vandermon, alag la cac he so Lagrange */ double p(double x); double p1(double x); void xemsl(double *x,double *y,int n); void nsvandermon(double *a); /*xqs[i] la cac diem quan sat, yqs[i] la cac ket qua quan sat, a[0],a[1],...,a[nqs] la cac he so can xac dinh cua da thuc giai truc tiep bang ma tran Vandermon*/ void noisuyham(double (*f)(double),double (*f1)(double),double a,double b); //Noi suy ham f tren khoang [a,b] void noisuyhamdb(); double poli(double x); int gjordan(kmatran a,double *x, int n); /*Dua ma tran a ve dang ma tran don vi. Giai he phuong trinh tuyen tinh bang khu gauss, tra ve true neu co nghiem */ //=================================================================== /*Tra ve gia tri da thuc noi suy tai x; avan[i] la cac he so cua da thuc giai truc tiep tu ma tran Vandermon, xqs[i] la cac diem quan sat*/ double poli(double x) //Tinh da thuc bang phuong phap Horner {int i;double s; s=avan[nqs]; for(i=nqs-1;i>=0;i--) s= s*x+avan[i]; return s; } //=============================================== /*Tra ve dao ham gia tri da thuc noi suy tai x; avan[i] la cac he so cua da thuc giai truc tiep tu ma tran Vandermon, xqs[i] la cac diem quan sat*/ double poli1(double x) //Tinh da thuc bang phuong phap Horner {int i;double s; s=nqs*avan[nqs]; for(i=nqs-1;i>0;i--) s= s*x+i*avan[i]; return s; } //=============================================== void xemsl(double *x,double *y,int n) {clrscr(); if(nqs<2) {cout<<"\nCo le ban chua nhap so lieu!";delay(2000);return;} int i; gotoxy(1,2); cout<<"TINH XAP XI DAO HAM"; cout<<endl<<"So mau quan sat n = "<<n; cout<<endl<<" i x[i] y[i]"<<endl; for(i=0;i<=n;i++) {cout<<endl<<setw(2)<<i; cout<<setw(8)<<setprecision(2)<<x[i]; cout<<setw(8)<<setprecision(2)<<y[i]; } gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); } //=============================================== /*Noi suy bang cach giai truc tiep he phuong trinh tuyen tinh voi ma tran Vandermon */ void nsvandermon(double *a) {int i,j,k,n1;kmatran aa;kvecto b; //Tinh ma tran Vandermon for(i=0;i<=nqs;i++) {aa[i][0]=1; for(j=1;j<=nqs;j++) aa[i][j]=aa[i][j-1]*xqs[i]; } for(i=0;i<=nqs;i++) aa[i][nqs+1]=yqs[i]; gjordan(aa,a,nqs); } //=============================================== void noisuyham(double (*f)(double),double (*f1)(double),double a,double b) {int i,j,nn,ntmp;double h;kvecto xx,yy;char ordst[2]; /*Gan cac gia tri xqs[i],yqs[i] va nqs cho xx[i],yy[i] va nn de khoi phuc lai truoc khi thoat khoi ham nay*/ nn=nqs; for(i=0;i<=nqs;i++) {xx[i]=xqs[i];yy[i]=yqs[i];} while(true) {clrscr(); cout>ntmp; if(ntmp<1) {cout<<"Bac da thuc phai khong am";delay(1000);continue;} if(ntmp>30) {cout<<"Bac da thuc qua lon";delay(1000);continue;} clrscr(); break; } nqs=ntmp; char *kieuham; if(f==sin) kieuham="Noi suy tren ham sin(x)";else kieuham=" "; //Ve do thi double scalex,scaley; /*[a,b] la khoang chua cac gia tri x, scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx,gocy,n,color; /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dijnh tuy thuoc vao dang ham*/ int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); while(true) {h=(b-a)/nqs; for(i=0;i<=nqs;i++) {xqs[i]=a+i*h;yqs[i]=f(xqs[i]);} nsvandermon(avan); for(i=nqs;i>=0;i--) if(avan[i]!=0) break; //Xac dinh bac cua da thuc noi suy Newton itoa(i,ordst,10);//Chuyen sang dang ky tu char ftitle[]="Da thuc noi suy Newton bac "; strcat(ftitle,ordst);//ftitle se la dong chu co them bac cua da thuc cleardevice(); setbkcolor(GREEN); setcolor(BLUE); gocx=gocy=0;//De chhuong trinh tu xac dinh n=200;//200 diem chia //Ve do thi ham f(x) vedothi(f,a,b,n,gocx,gocy,scalex,scaley,RED,1); //Ve do thi dao ham f'(x)=f1(x) vedothi(f1,a,b,n,gocx,gocy,scalex,scaley,MAGENTA,1); double aa=xqs[0],bb=xqs[nqs]; //Ve do thi da thuc noi suy vedothi(poli,aa,bb,n,gocx,gocy,scalex,scaley,WHITE,0); //Ve do thi dao ham da thuc noi suy vedothi(poli1,aa,bb,n,gocx,gocy,scalex,scaley,YELLOW,0); setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,kieuham); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve cac moc noi suy double x,y; for(i=0;i<=nqs;i++) {x=xqs[i]*scalex; y=-yqs[i]*scaley; setcolor(YELLOW); setfillstyle(SOLID_FILL,YELLOW); circle(x,0,2); floodfill(x,0,YELLOW); circle(x,y,2); floodfill(x,y,YELLOW); delay(5); } //Ve chu thich ham va da thuc noi suy settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); setcolor(RED); int x1,y1,x2,y2; x1=(getmaxx()-gocx)/2; y1=30-gocy; x2=x1+50; y2=y1; line(x1,y1,x2,y2); line(x1,y1+1,x2,y2+1);//Ve cho dam hon outtextxy(x2+10,y2,"Ham"); y1=y1+20;y2=y2+20; setcolor(WHITE); line(x1,y1,x2,y2); //line(x1,y1+1,x2,y2+1); outtextxy(x2+10,y2,"Da thuc ns"); y1=y1+20;y2=y2+20; setcolor(MAGENTA); line(x1,y1,x2,y2); line(x1,y1+1,x2,y2+1); outtextxy(x2+10,y2,"Dao ham"); y1=y1+20;y2=y2+20; setcolor(YELLOW); line(x1,y1,x2,y2); //line(x1,y1+1,x2,y2+1); outtextxy(x2+10,y2,"Dao ham dtns"); setcolor(BLUE); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"PgUp:"); outtextxy(x2+10,y2,"Tang 1 bac"); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"PgDn:"); outtextxy(x2+10,y2,"Giam 1 bac"); char mchon=getch(); if(mchon==0) mchon=getch(); if(mchon==73 && nqs<30)//Bam PgUp {nqs++;continue;} if(mchon==81 && nqs>0) //Bam PgDn {nqs--;continue;} break; } closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"\n\nTINH XAP XI DAO HAM DUNG DA THUC NOI SUY:"); fprintf(fp,"\n========================================"); fprintf(fp,"\n(%s)",kieuham); fprintf(fp,"\n(Da thuc noi suy bac %s)",ordst); fprintf(fp,"\n\nCAC DIEM MOC, HAM, DA THUC NOI SUY, DAO HAM, DHDTNS:"); fprintf(fp,"\n x[i] f(x[i]) dtns p(x) f'(x) p'(x)"); for(i=0;i<nqs;i++) fprintf(fp,"\n%10.4f %10.4f %10.4f %10.4f %10.4f",xqs[i],f(xqs[i]), poli(xqs[i]),f1(xqs[i]),poli1(xqs[i])); h=(xqs[1]-xqs[0])/2; fprintf(fp,"\n\nCAC DIEM GIUA MOC, HAM, DA THUC NOI SUY, DAO HAM, DHDTNS:"); fprintf(fp,"\n x[i]+h f(x[i]+h) dtns p(x) f'(x) p'(x)"); for(i=0;i<nqs;i++) fprintf(fp,"\n%10.4f %10.4f %10.4f %10.4f %10.4f",xqs[i]+h,f(xqs[i]+h), poli(xqs[i]+h),f1(xqs[i]+h),poli1(xqs[i]+h)); fclose(fp); xemtep(resultfile); gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); nqs=nn; for(i=0;i<=nqs;i++) {xqs[i]=xx[i];yqs[i]=yy[i];} //Tra lai gia tri cho nqs,xqs[i] va yqs[i] } //=============================================== void main() {int i,j;char mchon;double a,b; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(20,2);cout<<"TINH DAO HAM DUNG DA THUC NOI SUY"; gotoxy(2,2); cout<<endl<<" 1. Xem so lieu"; cout<<endl<<" 2. Dao ham sin(x)"; cout<<endl<<" 3. Dao ham 1/(1+x*x)"; cout<<endl<<" z. Ket thuc"; cout z de chon"; mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':xemsl(xqs,yqs,nqs); break; case '2':clrscr(); a=-2*M_PI;b=2*M_PI; noisuyham(sin,cos,a,b); break; case '3':clrscr(); a=0;b=1; noisuyham(p,p1,a,b); break; } //gotoxy(50,25); //cout<<"Nhan phim bat ky de tiep tuc"; //getch(); } } //=============================================== /*Giai he phuong trinh tuyen tinh bang cach dua ma tran a ve dang ma tran don vi. Tra ve gia tri true neu co nghiem */ int gjordan(kmatran a,double *x,int n) {int i,j,k,h;double tmp,p; int n1=n+1; for(i=0;i<=n;i++) //Vong lap cac buoc khu {//Tim hang co phan tu dau lon nhat h=i; for(k=i+1;k<=n;k++) if(fabs(a[k][i])> fabs(a[h][i])) h=k; if(a[h][i]==0) {cout<<"Ma tran suy bien";delay(1000);return false;} if(h!=i) //Doi hang i va hang h vi a[h][i] > a[i][i] {int j;double tmp; for(j=i;j<=n1;j++) {tmp=a[i][j];a[i][j]=a[h][j];a[h][j]=tmp;} } //chuyen he so a[i][i] = 1 tmp=a[i][i]; for(j=i;j<=n1;j++) a[i][j] = a[i][j]/tmp; //Bat tinh lai cac hang for(k=0;k<=n;k++) {if(k==i) continue; p=a[k][i]; /*Vi ta biet a[k][i] =0 sau bien doi, chi tinh tu a[k][i+1]*/ for(j=i+1;j<=n1;j++) a[k][j]=a[k][j] - p*a[i][j]; } } for(i=0;i<=n;i++) x[i]=a[i][n+1]; return true; } //=============================================== double p(double x) {return(1/(1+x*x)); } //=============================================== double p1(double x) {return(-2*x/((1+x*x)*(1+x*x))); } Phương pháp hình thang tính tích phân xác định: HTHANG.CPP //HTHANG.CPP /*Phuong phap tinh xap xi tich phan bang phuong phap hinh thang tren khoang [a,b]*/ #include "zheader.cpp" #include "zlibrary.cpp" char *ftitle,*datafile="HTHANG.DAT",*resultfile="HTHANG.KQ"; typedef double kvecto[200]; double epsi=0.01; int nmax=100; double dlc=1,m=0; //=============================================== double pi(double x) {return (4/(1+x*x)); } //=============================================== //a=-2;b=1; double gauss(double x) {return (exp(-(x-m)*(x-m)/(2*dlc*dlc)))/(dlc*sqrt(2*M_PI)); } //=============================================== int hinhthang(double (*f)(double),double a,double b,double >tp, double &err,int &khoangchia) /*Phuong phap hinh thang tinh tich phan xac dinh trong khoang [a,b]. Bien gttp la gia tri xap xi cua tich phan tinh duoc. Tra ve gia tri true neu da dat duoc do chinh xac*/ {clrscr(); {int c=2,h=2; gotoxy(c,h); printf("CAC THONG SO DE TINH TOAN:"); gotoxy(c,++h); printf("Can duoi a = : %6.1f",a); gotoxy(c,++h); printf("Can tren b = : %6.1f",b); gotoxy(c,++h); printf("Khoang chia toi da: %2d",nmax); gotoxy(c,++h); printf("Sai so : %6.5f",epsi); gotoxy(c,h+2); printf("Co chap nhan khong (C/K)"); char mchon=toupper(getch()); if(mchon=='K') {h+=2; gotoxy(c,h); printf("NHAP LAI CAC THONG SO: "); gotoxy(c,++h); printf("Can duoi a = :");scanf("%lf",&a); gotoxy(c,++h); printf("Can tren b = :");scanf("%lf",&b); gotoxy(c,++h); printf("Khoang chia toi da:");scanf("%d",&nmax); gotoxy(c,++h); printf("Sai so :");scanf("%lf",&epsi); } } clrscr(); double s0,s1,s2,h,tp,tp1;int k,nkc,i; kvecto x; nkc=1; h=b-a; s0=(f(a)+f(b))/2; s1=0; tp=(s0+s1)*h; do {tp1=tp; nkc=nkc*2; h=h/2; s2=0;//Bat dau tinh tong tai cac diem moi for(i=1;i<nkc;i+=2) s2+=f(a+i*h); s1=s1+s2; tp=h*(s0+s1); if(nkc>nmax) {cout<<endl<<"Tich phan chua hoi tu voi "<<nmax<<" khoang chia"; delay(1000);return false; } } while(fabs(tp-tp1)>epsi); err=fabs(tp-tp1);khoangchia=nkc;gttp=tp; for(i=0;i<=nkc;i++) x[i]=a+h*i; //Ve do thi double scalex,scaley; /*scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx=0,gocy=0,n=100,color=BLUE;//100 diem chia /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dinh tuy thuoc vao dang ham*/ int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); cleardevice(); setbkcolor(LIGHTGRAY); setcolor(BLUE); vedothi(f,a,b,n,gocx,gocy,scalex,scaley,BLUE,1); //Ve cac diem int x1,y1,x2,y2,d; color=WHITE; setcolor(color); setfillstyle(SOLID_FILL,YELLOW); for(i=0;i<nkc;i++) {x1=x[i]*scalex; y1=-f(x[i])*scaley; x2=x[i+1]*scalex; y2=-f(x[i+1])*scaley; line(x1,0,x1,y1); line(x2,0,x2,y2); line(x1,y1,x2,y2); circle(x1,0,2); floodfill(x1,0,color); delay(5); } circle(x2,0,2); floodfill(x2,0,color); //Ve cac thong bao ket qua setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,"PHUONG PHAP HINH THANG"); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve chu thich ham settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); char st1[20],st2[10],st3[10],st4[10],st5[10]; sprintf(st1,"[%4.1f,%4.1f]",a,b); sprintf(st2,"%6.4f",gttp); sprintf(st3,"%2d",khoangchia); sprintf(st4,"%6.4f",err); x1=(getmaxx()-gocx)/2-25; y1=30-gocy; x2=x1+65; y2=y1; outtextxy(x1,y1,"[a,b] :"); outtextxy(x2+10,y2,st1); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"T/phan :"); outtextxy(x2+10,y2,st2); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"K/chia :"); outtextxy(x2+10,y2,st3); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Sai so :"); outtextxy(x2+10,y2,st4); getch(); closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"XAP XI TICH PHAN XAC DINH BANG P/P HINH THANG"); fprintf(fp,"\n(%s)",ftitle); fprintf(fp,"\n\nKET QUA TINH TOAN:"); fprintf(fp,"\n==================="); fprintf(fp,"\nTich phan tren khoang [%4.1f,%4.1f]",a,b); fprintf(fp,"\nGia tri xap xi tich phan: %6.4f",gttp); fprintf(fp,"\nSo khoang chia: %d",khoangchia); fprintf(fp,"\nSai so giua 2 gia tri xap xi cuoi cung: %6.4f",err); fclose(fp); xemtep(resultfile); return true; } //=============================================== void main() {int kchia;double a,b,tpxd,err; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(2,2); cout<<"XAP XI TICH PHAN XAC DINH BANG PHUONG PHAP HINH THANG"; cout<<endl; cout<<endl<<" 1. sin(x)";//a=-2,b=2; cout<<endl<<" 2. exp(x)";//a=0,b=1; cout<<endl<<" 3. Mat do phan bo chuan";//a=-2,b=1; cout<<endl<<" 4. 4/(1+x*x) (so pi)";//a=0,b=1; cout<<endl<<" z. Khong tinh nua"; cout z de chon"; char mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr(); a=-2;b=3; ftitle="Xap xi tich phan ham sin(x)"; hinhthang(sin,a,b,tpxd,err,kchia); break; case '2':clrscr(); a=0;b=1; ftitle="Xap xi tich phan ham exp(x)"; hinhthang(exp,a,b,tpxd,err,kchia); break; case '3':clrscr(); a=-2;b=1; ftitle="Xap xi tich phan ham mat do phan bo chuan"; hinhthang(gauss,a,b,tpxd,err,kchia); break; case '4':clrscr(); a=0;b=1; ftitle="Xap xi tich phan ham 4/(1+x*x) (=so pi)"; hinhthang(pi,a,b,tpxd,err,kchia); break; } gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); } } Tính tích phân bằng phương pháp Simson: SIMSON.CPP //SIMSON.CPP /*Phuong phap tinh xap xi tich phan bang phuong phap Simson (pp cau phuong) tren khoang [a,b]*/ #include "zheader.cpp" #include "zlibrary.cpp" char *ftitle,*datafile="SIMSON.DAT",*resultfile="SIMSON.KQ"; typedef double kvecto[200]; double epsi=0.5; int nmax=100; double dlc=1,m=0; double a0,a1,a2,xqs0,xqs1;//De tinh da thuc bac 2 //============================================= double pi(double x) {return (4/(1+x*x)); } //============================================= //a=-2;b=1; double gauss(double x) {return (exp(-(x-m)*(x-m)/(2*dlc*dlc)))/(dlc*sqrt(2*M_PI)); } //============================================= /*Tinh lai cac he so da thuc bac 2 a0,a1,a2 di qua 3 diem (x0,f(x0)), (x1,f(x1)),(x2,f(x2))*/ void nsnewton(double x0,double y0,double x1,double y1,double x2,double y2) {double sp0,sp1; a0=y0; sp0=(y1-y0)/(x1-x0); sp1=(y2-y1)/(x2-x1); a1=sp0; a2=(sp1-sp0)/(x2-x0); } //============================================= double newpoli(double x) {return(a0+a1*(x-xqs0)+a2*(x-xqs0)*(x-xqs1)); } //============================================= int simson(double (*f)(double),double a,double b,double >tp, double &err,int &khoangchia) /*Phuong phap hinh thang tinh tich phan xac dinh trong khoang [a,b]. Bien gttp la gia tri xap xi cua tich phan tinh duoc. Tra ve gia tri true neu da dat duoc do chinh xac*/ {clrscr(); {int c=2,h=2; gotoxy(c,h); printf("CAC THONG SO DE TINH TOAN:"); gotoxy(c,++h); printf("Can duoi a = : %6.1f",a); gotoxy(c,++h); printf("Can tren b = : %6.1f",b); gotoxy(c,++h); printf("Khoang chia toi da: %2d",nmax); gotoxy(c,++h); printf("Sai so : %6.5f",epsi); gotoxy(c,h+2); printf("Co chap nhan khong (C/K)"); char mchon=toupper(getch()); if(mchon=='K') {h+=2; gotoxy(c,h); printf("NHAP LAI CAC THONG SO: "); gotoxy(c,++h); printf("Can duoi a = :");scanf("%lf",&a); gotoxy(c,++h); printf("Can tren b = :");scanf("%lf",&b); gotoxy(c,++h); printf("Khoang chia toi da:");scanf("%d",&nmax); gotoxy(c,++h); printf("Sai so :");scanf("%lf",&epsi); } } clrscr(); double s0,s1,s2,h,tp,tp1;int k,nkc,i; kvecto x; nkc=1; h=b-a; h=b-a; s0=f(a)+f(b); s1=0; s2=0; tp=(s0+2*s1+4*s2)*h/3; do {tp1=tp; s1=s1+s2; nkc=nkc*2; h=h/2; s2=0;//Bat dau tinh tong tai cac diem moi for(i=1;i<nkc;i+=2) s2=s2+f(a+i*h); tp=(s0+2*s1+4*s2)*h/3; if(nkc>nmax) {cout<<endl<<"Tich phan chua hoi tu voi "<<nmax<<" khoang chia"; delay(1000);return false; } } while(fabs(tp-tp1)>epsi); err=fabs(tp-tp1);khoangchia=nkc;gttp=tp; for(i=0;i<=nkc;i++) x[i]=a+h*i; //Ve do thi double scalex,scaley; /*scalex,scaley la do tang (hoac giam neu <1) cac gia tri theo truc hoanh va truc tung*/ int gocx=0,gocy=0,n=100,color=RED;//100 diem chia /*(gocx,gocy) la diem toa do. Neu toa do khac (0,0) thi cac gia tri toa do va scalex,scaley se giu nguyen, neu bang 0 thi chuong trinh se xac dinh tuy thuoc vao dang ham*/ int mh=DETECT,gmode=0; initgraph(&mh,&gmode,"\\TC\\BGI"); cleardevice(); setbkcolor(GREEN); setcolor(BLUE); vedothi(f,a,b,n,gocx,gocy,scalex,scaley,RED,1); //Ve cac diem int x1,y1,x2,y2; double aa,bb; color=WHITE; setcolor(color); setfillstyle(SOLID_FILL,YELLOW); for(i=0;i<nkc;i++) {x1=x[i]*scalex; y1=-f(x[i])*scaley; x2=x[i+1]*scalex; y2=-f(x[i+1])*scaley; line(x1,0,x1,y1); line(x2,0,x2,y2); circle(x1,0,2); floodfill(x1,0,color); delay(5); } circle(x2,0,2); floodfill(x2,0,color); //Ve cac da thuc bac 2 double xx0,yy0,xx1,yy1,xx2,yy2,nn=10; for(i=0;i<nkc-1;i+=2) {xx0=x[i]; yy0=f(x[i]); xx1=x[i+1]; yy1=f(x[i+1]); xx2=x[i+2]; yy2=f(x[i+2]); //Tinh lai cac he so da thuc nsnewton(xx0,yy0,xx1,yy1,xx2,yy2); aa=xx0;bb=xx2;xqs0=xx0;xqs1=xx1; vedothi(newpoli,aa,bb,nn,gocx,gocy,scalex,scaley,color,0); delay(5); } //Ve cac thong bao ket qua setcolor(BLUE); settextstyle(TRIPLEX_FONT,HORIZ_DIR,1); outtextxy(0,getmaxy()-gocy-40,"PHUONG PHAP SIMSON (CAU PHUONG)"); outtextxy(0,getmaxy()-gocy-20,ftitle); //Ve chu thich ham settextstyle(DEFAULT_FONT,HORIZ_DIR,1); settextjustify(LEFT_TEXT,CENTER_TEXT); char st1[20],st2[10],st3[10],st4[10],st5[10]; sprintf(st1,"[%4.1f,%4.1f]",a,b); sprintf(st2,"%6.4f",gttp); sprintf(st3,"%2d",khoangchia); sprintf(st4,"%6.4f",err); x1=(getmaxx()-gocx)/2-25; y1=30-gocy; x2=x1+65; y2=y1; outtextxy(x1,y1,"[a,b] :"); outtextxy(x2+10,y2,st1); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"T/phan :"); outtextxy(x2+10,y2,st2); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"K/chia :"); outtextxy(x2+10,y2,st3); y1=y1+20;y2=y2+20; outtextxy(x1,y1,"Sai so :"); outtextxy(x2+10,y2,st4); getch(); closegraph(); //Dua ket qua vao tep ketqua FILE *fp=fopen(resultfile,"wt"); fprintf(fp,"XAP XI TICH PHAN XAC DINH BANG P/P SIMSON"); fprintf(fp,"\n(%s)",ftitle); fprintf(fp,"\n\nKET QUA TINH TOAN:"); fprintf(fp,"\n==================="); fprintf(fp,"\nTich phan tren khoang [%4.1f,%4.1f]",a,b); fprintf(fp,"\nGia tri xap xi tich phan: %6.4f",gttp); fprintf(fp,"\nSo khoang chia: %d",khoangchia); fprintf(fp,"\nSai so giua 2 gia tri xap xi cuoi cung: %6.4f",err); fclose(fp); xemtep(resultfile); return true; } //============================================= void main() {int kchia;double a,b,tpxd,err; while(true) {textbackground(BLACK); textcolor(WHITE); clrscr(); gotoxy(2,2); cout<<"XAP XI TICH PHAN XAC DINH BANG PHUONG PHAP SIMSON"; cout<<endl; cout<<endl<<" 1. sin(x)";//a=-2,b=2; cout<<endl<<" 2. exp(x)";//a=0,b=1; cout<<endl<<" 3. Mat do phan bo chuan";//a=-2,b=1; cout<<endl<<" 4. 4/(1+x*x) (so pi)";//a=0,b=1; cout<<endl<<" z. Khong tinh nua"; cout z de chon"; char mchon=toupper(getch()); if(mchon=='Z') break; switch(mchon) {case '1':clrscr(); a=-2;b=3; ftitle="Xap xi tich phan ham sin(x)"; simson(sin,a,b,tpxd,err,kchia); break; case '2':clrscr(); a=0;b=1; ftitle="Xap xi tich phan ham exp(x)"; simson(exp,a,b,tpxd,err,kchia); break; case '3':clrscr(); a=-2;b=1; ftitle="Xap xi tich phan ham mat do phan bo chuan"; simson(gauss,a,b,tpxd,err,kchia); break; case '4':clrscr(); a=0;b=1; ftitle="Xap xi tich phan ham 4/(1+x*x) (=so pi)"; simson(pi,a,b,tpxd,err,kchia); break; } gotoxy(50,25); cout<<"Nhan phim bat ky de tiep tuc"; getch(); } } Tài liệu tham khảo 1. Dương Thùy Vỹ, Phương pháp tính, Nhà xuất bản Khoa học và Kỹ thuật, 2001 2. Đinh Văn Phong, Phương pháp số trong cơ học, Nhà xuất bản Khoa học và Kỹ thuật, 1999 3. Lê Trọng Vinh, Giải tích số, Nhà xuất bản Khoa học và Kỹ thuật, 2000 4. Phạm Kỳ Anh , Giải tích số, Nhà xuất bản Đại học Quốc gia Hà Nội, 1966 5. Phạm Phú Triêm - Nguyễn Bường, Giải tích số, thuật toán, chương trình Pascal, Nhà xuất bản Đại học Quốc gia Hà Nội, 2000 6. Szidarovszky Ferenc, Nhập môn phương pháp số (tiếng Hung), Kozgazdasági és Jogi Konyvkiado, Budapest 1974. 7. Tạ Văn Đỉnh, Phương pháp tính, Nhà xuất bản Giáo dục - 1995 8. Trần Văn Minh, Phương pháp số và chương trình bằng Turbo Pascal, Nhà xuất bản Khoa học và Kỹ thuật, 1998 9. X.M Ermakov, Phương pháp Monte-Carlo và các vấn đề liên quan, Nhà xuất bản Khoa học và Kỹ thuật, 1977