#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
