编辑代码

#include <stdio.h>
int main ()
 {
   double x;
   double x1,x2;
   double xmin,xmax;
   double k1,k2,m1,m2,p1,p2;//k=a-x,m=b-nx,p=m^n≈m^2,这里约取n=2
   double a1,b1,a2,b2;//这四个符号对应书上的a1,b1,a2,b2
   double y,y0,y1,y2;//y1,y2=k*p=(a-x)*(b-2*x)^2;y,y0=y1-y2
   int s=0;
   a1=8.9289; b1=1.0711; a2=3.8605; b2=1.1395;//默认a1<a2的情况
   x1=0; x2=b1;//x1设为0，x2设为b1,作为初始区间(0,b1),另外取x2=b2不影响结果
   //这里取x2为b1是因为a比b大，实际上结果得到的x是公式里的“2x”
   //如果b比a大，应该取x2=a1或a2，并把下面的0.5*xmin、xmax前面的0。5*去掉，在后面的m=b-x改成m=b-2*x;
   while(s<100)//循环100次，但实际上大概20次之前结果就已经出来了
   {
    s=s+1;
    x=(x1+x2)/2;
    xmin=(x1+x)/2;xmax=(x2+x)/2;
    k1=a1-0.5*xmin;m1=b1-xmin;p1=m1*m1;y1=k1*p1;//取公式中的n为2
    k2=a2-0.5*xmin;m2=b2-xmin;p2=m2*m2;y2=k2*p2;
    y=y1-y2;
    k1=a1-0.5*xmax;m1=b1-xmax;p1=m1*m1;y1=k1*p1;
    k2=a2-0.5*xmax;m2=b2-xmax;p2=m2*m2;y2=k2*p2;
    y0=y1-y2;
   printf("%f %f %f %f %f %f %f %d\n",y1,y2,y,y0,x,x1,x2,s);
   if(y*y>y0*y0) x1=x;
   if(y*y<y0*y0) x2=x;
   }
   //试求y=(a1-x)*(b1-nx)^n-(a2-x)*(b2-nx)^n=0在(0,a1)上的解
   //对于区间(x1,x2)，取x为x1与x2平均值，再分别取xmin、xmax为x与x1、x2的平均值
   //分别求xmin与xmax对应的y=(a1-x)*(b1-nx)^n-(a2-x)*(b2-nx)的绝对值更小
   //并把对应的x换成x1或者x2
   //例如，如果用xmin算出来的y绝对值更小，就把x换成x2，(x1,x2)→(x1,x)(即(x1,1/2(x1+x2))
   //因为xmin算出来的y绝对值更小说明解x0应该落在(x1,x)这半边
   printf("%f",x);
	return 0;
}