DFP.rar_C Builder_dfp

#include<stdio.h> #include<conio.h> #include<math.h> #include<stdlib.h> float H[2][2]={{1,0},{0,1}},HT[2][2], grad[1][2], gradT[2][1], s[2][1],sT[1][2],g[2][1],gT[1][2], rezpro[2][2]; float delta=0.001; int menu(); void mpro(float a[][2], float b[][2], int n, int m, int t); void mpro(float a[][2], float b[][1], int n, int m, int t); void mpro(float a[2][1], float b[1][2], int n, int m, int t); void napr(); float planer(float x10, float x20); float sven(float x, float x1, float x2); float roz(float c, float d); float rozplane(float t, float x1, float x2); void gradRoz(float a, float b); void peres(float x1, float x2, int id); float norma(); void gradKvad(float a, float b); float kva(float c, float d); float kvaplane(float t, float x1, float x2); float svenk(float x, float x1, float x2); float planek(float x10, float x20); void main(void) { float x10, x20, tm; int w, o=0; FILE* ro; FILE* kv; a1: H[0][0]=1; H[1][0]=0; H[0][1]=0; H[1][1]=1; w=menu(); o=0; switch(w) { case 1: ro=fopen("Roz.txt","w"); printf("X10="); scanf("%f",&x10); printf("X20="); scanf("%f",&x20); gradRoz(x10,x20); while(norma()>0.01){ o++; napr(); tm=planer(x10,x20); x10+=s[0][0]*tm; x20+=s[1][0]*tm; s[0][0]=s[0][0]*tm; s[1][0]=s[1][0]*tm; peres(x10,x20, 0); printf("iteration %d\n",o); printf("X1=%f\n",x10); printf("X2=%f\n\n",x20); fprintf(ro,"iteration %d\n",o); fprintf(ro,"X1=%f\n",x10); fprintf(ro,"X2=%f\n\n",x20); } printf("Norma gradienta: %f\n",norma()); fprintf(ro,"Norma gradienta: %f",norma()); getch(); fclose(ro); break; case 2: kv=fopen("Kvad.txt","w"); printf("X10="); scanf("%f",&x10); printf("X20="); scanf("%f",&x20); gradKvad(x10,x20); while(norma()>0.01){ o++; napr(); tm=planek(x10,x20); x10+=s[0][0]*tm; x20+=s[1][0]*tm; s[0][0]=s[0][0]*tm; s[1][0]=s[1][0]*tm; peres(x10,x20, 1); printf("iteration %d\n",o); printf("X1=%f\n",x10); printf("X2=%f\n\n",x20); fprintf(kv,"iteration %d\n",o); fprintf(kv,"X1=%f\n",x10); fprintf(kv,"X2=%f\n\n",x20); } printf("Norma gradienta: %f\n",norma()); fprintf(kv,"Norma gradienta: %f\n",norma()); getch(); break; case 3: exit(1); } goto a1; } int menu() { int q; printf("choose function type\n\n1-Rozenbrok y=0.1*(x2-x1^2)^2+(5-x1)^2\n\n2-kvadratichnaya y=0.1*x1^2+0.5*x1*x2+5*x2^2\n\n3-exit\n\n"); scanf("%d",&q); return q; } void gradRoz(float a, float b) { grad[0][0]=-0.1*4*a*(b-a*a)-2*(5-a); grad[0][1]=0.1*2*(b-a*a); } float roz(float c, float d) { float y; y=0.1*(d-c*c)*(d-c*c)+(5-c)*(5-c); return y; } float rozplane(float t, float x1, float x2) { float y; y=0.1*((x2+s[1][0]*t)-(x1+s[0][0]*t)*(x1+s[0][0]*t))*((x2+s[1][0]*t)-(x1+s[0][0]*t)*(x1+s[0][0]*t))+(5-(x1+s[0][0]*t))*(5-(x1+s[0][0]*t)); return y; } void mpro(float a[][2], float b[][2], int n, int m, int t) { float q; int i, j, k; for(k=0;k<n;k++) { for(j=0;j<t;j++) { q=0; for(i=0;i<m;i++) { q=q+a[k][i]*b[i][j]; } rezpro[k][j]=q; } } } void mpro(float a[][2], float b[][1], int n, int m, int t) { float q; int i, j, k; for(k=0;k<n;k++) { for(j=0;j<t;j++) { q=0; for(i=0;i<m;i++) { q=q+a[k][i]*b[i][j]; } rezpro[k][j]=q; } } } void mpro(float a[2][1], float b[1][2], int n, int m, int t) { float q; int i, j, k; for(k=0;k<n;k++) { for(j=0;j<t;j++) { q=0; for(i=0;i<m;i++) { q=q+a[k][i]*b[i][j]; } rezpro[k][j]=q; } } } void napr() { int i; for(i=0;i<2;i++) { gradT[i][0]=grad[0][i]; } mpro(H,gradT,2,2,1); for(i=0;i<2;i++) { s[i][0]=-rezpro[i][0]; } } float planer(float x10, float x20) { float t1=x10, t2, tmin, pr=1; t1=0; t2=sven(0, x10, x20); if(t1<t2) { pr=(rozplane(t1+delta,x10,x20)-rozplane(t1-delta,x10,x20))/(2*delta); tmin=-0.5*pr/(rozplane(t2, x10, x20)-rozplane(t1, x10, x20)-pr); } else { pr=(rozplane(t2+delta,x10,x20)-rozplane(t2-delta,x10,x20))/(2*delta); tmin=-0.5*pr/(rozplane(t1, x10, x20)-rozplane(t2, x10, x20)-pr); } return tmin; } float sven(float x, float x1, float x2) { float pr1, pr2; float del=1, xx; if(rozplane(x, x1, x2)<rozplane(x+del,x1,x2)) { pr1=(rozplane(x+delta,x1,x2)-rozplane(x-delta,x1,x2))/(2*delta); pr2=(rozplane(x+del+delta,x1,x2)-rozplane(x+del-delta,x1,x2))/(2*delta); if(pr1*pr2>0) { del*=-1; pr2=(rozplane(x+del+delta,x1,x2)-rozplane(x+del-delta,x1,x2))/(2*delta); if(pr1*pr2<0) { return x+del; } } else { return x+del; } } xx=x+del; pr1=(rozplane(x+delta,x1,x2)-rozplane(x-delta,x1,x2))/(2*delta); pr2=(rozplane(xx+delta,x1,x2)-rozplane(xx-delta,x1,x2))/(2*delta); while(pr1*pr2>0) { del*=2; xx+=del; pr2=(rozplane(xx+delta,x1,x2)-rozplane(xx-delta,x1,x2))/(2*delta); } return xx; } void peres(float x1, float x2, int id) { int i, j; float a[2][2], b, c[2][1], d[2][2], e[2][2], f[1][2], r; float opr; for(i=0;i<2;i++) { g[i][0]=grad[0][i]; sT[0][i]=s[i][0]; } switch(id) { case 0: gradRoz(x1,x2); break; case 1: gradKvad(x1,x2); break; } g[0][0]=grad[0][0]-g[0][0]; g[1][0]=grad[0][1]-g[1][0]; for(i=0;i<2;i++) { gT[0][i]=g[i][0]; } opr=H[0][0]*H[1][1]-H[0][1]*H[1][0]; HT[0][0]=H[1][1]/opr; HT[1][1]=H[0][0]/opr; HT[0][1]=-H[1][0]/opr; HT[1][0]=-H[0][1]/opr; mpro(s,sT,2,1,2); for(i=0;i<2;i++) { for(j=0;j<2;j++) { a[i][j]=rezpro[i][j]; rezpro[i][j]=0; } } mpro(sT,g,1,2,1); b=rezpro[0][0]; rezpro[0][0]=0; mpro(H,g,2,2,1); for(i=0;i<2;i++) { c[i][0]=rezpro[i][0]; rezpro[i][0]=0; } mpro(c,gT,2,1,2); for(i=0;i<2;i++) { for(j=0;j<2;j++) { d[i][j]=rezpro[i][j]; rezpro[i][j]=0; } } mpro(d,H,2,2,2); for(i=0;i<2;i++) { for(j=0;j<2;j++) { e[i][j]=rezpro[i][j]; rezpro[i][j]=0; } } mpro(gT,HT,1,2,2); for(i=0;i<2;i++) { f[0][i]=rezpro[0][i]; rezpro[0][i]=0; } mpro(f,g,1,2,1); r=rezpro[0][0]; rezpro[0][0]=0; for(i=0;i<2;i++) { for(j=0;j<2;j++) { H[i][j]=H[i][j]+a[i][j]/b-e[i][j]/r; } } } float norma() { float q; q=sqrt(grad[0][0]*grad[0][0]+grad[0][1]*grad[0][1]); return q; } void gradKvad(float a, float b) { grad[0][0]=0.1*2*a+0.1*5*b; grad[0][1]=5*2*b+0.1*5*a; } float kva(float c, float d) { float y; y=0.1*c*c+0.1*5*c*d+5*d*d; return y; } float kvaplane(float t, float x1, float x2) { float y; y=0.1*((x1+s[0][0]*t)*(x1+s[0][0]*t))+0.1*5*((x1+s[0][0]*t)*(x2+s[1][0]*t))+5*((x2+s[1][0]*t)*(x2+s[1][0]*t)); return y; } float planek(float x10, float x20) { float t1=x10, t2, tmin, pr=1; t1=1; t2=svenk(t1, x10, x20); if(t1<t2) { pr=(kvaplane(t1+delta,x10,x20)-kvaplane(t1-delta,x10,x20))/(2*delta); tmin=-0.5*pr/(kvaplane(t2, x10, x20)-kvaplane(t1, x10, x20)-pr); } else { pr=(kvaplane(t2+delta,x10,x20)-kvaplane(t2-delta,x10,x20))/(2*delta); tmin=-0.5*pr/(kvaplane(t1, x10, x20)-kvaplane(t2, x10, x20)-pr); } return tmin; } float svenk(float x, float x1, float x2) { float pr1, pr2; float del=1, xx; if(kvaplane(x, x1, x2)<kvaplane(x+del,x1,x2)) { pr1=(kvaplane(x+delta,x1,x2)-kvaplane(x-delta,x1,x2))/(2*delta); pr2=(kvaplane(x+del+delta,x1,x2)-kvaplane(x+del-delta,x1,x2))/(2*delta); if(pr1*pr2>0) { del*=-1; pr2=(kvaplane(x+del+delta,x1,x2)-kvaplane(x+del-delta,x1,x2))/(2*delta); if(pr1*pr2<0) { return x+del; } } else { return x+del; } } xx=x+del; pr1=(kvaplane(x+delta,x1,x2)-kvaplane(x-delta,x1,x2))/(2*delta); pr2=(kvaplane(xx+delta,x1,x2)-kvaplane(xx-delta,x1,x2))/(2*delta); while(pr1*pr2>0) { del*=2; xx+=del; pr2