/*最小二乘法用抛物线拟合五个点的数值分析问题2006-11-17 21:49TC环境下编译通过.大家可以自己设置路径.题目如下:用一条抛物线y=ax*x+b*x+c拟合以下数据,使误差平方和最小        i     1     2     3     4     5
                           xi    0     4     8    24    40
                           yi    80   130  180     335    377
求出真解,

再画出 y=80+300sin(圆周率*x/74)图形.)以上五个点均在此曲线上)
与上面的拟合曲线画在同一坐标系下.标出真解*/

#include <stdio.h>
#include <graphics.h>
#include <conio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415927
double x[5]={0, 4 ,8, 24, 40};
double y[5]={80, 130, 180, 335, 377};

void draw1()
{
     int i, j;
     double x, y;

     setcolor(1);
     for (i=0; i<40; i++)
     {
          x=(double)i;
          y=80+300*sin(PI/74.0*x);
          if (i==0)
           moveto(10*x+100, getmaxy()-y);
          else
           lineto(10*x+100, getmaxy()-y);
     }
}

void draw2()
{
     int i, j;
     double p1[4]={0}, p2[4]={0}, p3[4]={0};
     double a, b, c;
     double x1, y1,p20,p30,p31;

     for (i=0; i<5; i++)
     {
          p1[0]+=pow(x[i], 4);
          p1[1]+=pow(x[i], 3);
          p1[2]+=pow(x[i], 2);
          p1[3]+=pow(x[i], 2)*y[i];
     }

     for (i=0; i<5; i++)
     {
          p2[0]+=pow(x[i], 3);
          p2[1]+=pow(x[i], 2);
          p2[2]+=x[i];
          p2[3]+=x[i]*y[i];
     }

     for (i=0; i<5; i++)
     {
          p3[0]+=pow(x[i], 2);
          p3[1]+=x[i];
          p3[2]+=1;
          p3[3]+=y[i];
     }

     p20=p2[0];
     p30=p3[0];
     for (i=0; i<4; i++)
     {
          p2[i]-=(p1[i]*p20/p1[0]);
          p3[i]-=(p1[i]*p30/p1[0]);
     }

     p31=p3[1];
     for (i=0; i<4; i++)
     {
            p3[i]-=(p2[i]*p31/p2[1]);
     }

     c=p3[3]/p3[2];
     b=(p2[3]-p2[2]*c)/p2[1];
     a=(p1[3]-p1[2]*c-p1[1]*b)/p1[0];

     printf("\na=%lf, b=%lf, c=%lf\n", a, b, c);

     setcolor(2);
     for (i=0; i<40; i++)
     {
          x1=(double)i;
          y1=a*x1*x1+b*x1+c;
          if (i==0)
           moveto(10*x1+100, getmaxy()-y1);
          else
           lineto(10*x1+100, getmaxy()-y1);
     }
}

void zuobiao()
{
     setcolor(15);
     line(100, getmaxy(), 550, getmaxy());
     line(140, getmaxy(), 140, getmaxy()-10);
     outtextxy(140, getmaxy()-20, "4");
     line(180, getmaxy(), 180, getmaxy()-10);
     outtextxy(180, getmaxy()-20, "8");
     line(340, getmaxy(), 340, getmaxy()-10);
     outtextxy(340, getmaxy()-20, "24");
     line(500, getmaxy(), 500, getmaxy()-10);
     outtextxy(500, getmaxy()-20, "40");

     line(100, getmaxy(), 100, 100);
     line(100, getmaxy()-80, 90, getmaxy()-80);
     outtextxy(65, getmaxy()-80, "80");
     line(100, getmaxy()-130, 90, getmaxy()-130);
     outtextxy(65, getmaxy()-130, "130");
     line(100, getmaxy()-180, 90, getmaxy()-180);
     outtextxy(65, getmaxy()-180, "180");
     line(100, getmaxy()-335, 90, getmaxy()-335);
     outtextxy(65, getmaxy()-335, "335");
     line(100, getmaxy()-377, 90, getmaxy()-377);
     outtextxy(65, getmaxy()-377, "377");

     outtextxy(300, 250, "BLUE: y=80+300*sin(PI*x/74)");
     outtextxy(390, 275, "2");
     outtextxy(300, 280, "GREEN: y=ax +bx+c");
}

void main()
{
     int gdriver=DETECT,gmode,errorcode;

     initgraph(&gdriver,&gmode,"D:\\TC\\bgi");
     errorcode=graphresult();

     if(errorcode!=grOk)
     {
         printf("Graphics error:%s\n",grapherrormsg(errorcode));
         printf("Press any key to halt:");
         getch();
         exit(1);
     }

     draw1();
     draw2();
     zuobiao();
 
     getch();
     closegraph();
}