用C语言实现的曲线拟合的最小二乘法(共10页).doc

精选优质文档-倾情为你奉上实验名称：曲线拟合的最小二乘法实验目的了解曲线拟合的最小二乘法实验类型设计型实验环境Windows XP TC实验内容相关知识：已知Ca,b中函数f(x)的一组实验数据（xi,yi）（i=0,1,m），其中yi=f(xi)。设是Ca,b上线性无关函数族。在中找函数f(x) 曲线拟合的最小二乘解，其法方程（组）为： 其中， k=0,1,n特别是，求函数f(x) 曲线拟合的线性最小二乘解的计算公式为：数据结构：两个一维数组或一个二维数组算法设计：（略）实验用例： 已知函数y=f(x)的一张表：x0102030405060708090y6867.166.465.664.661.861.060.860.460试验要求：利用曲线拟合的线性最小二乘法求被逼近函数f(x)在点x=55处的近似值，并画出实验数据和直线。编写代码：#include<stdio.h>#include<stdlib.h> #include<graphics.h>double qiuhe1(double a102,int p)int i;double y;y=0;for(i=0;i<10;i+)y=y+aip;return y;double qiuhe2(double a102,int p)int i;double y=0;for(i=0;i<10;i+)y=y+ai0*aip;return y;double nihe(double a102,double x)double a1,b,y;a1=(10*qiuhe2(a,1)-qiuhe1(a,0)*qiuhe1(a,1)/(10*qiuhe2(a,0)-qiuhe1(a,0)*qiuhe1(a,0);b=(qiuhe2(a,0)*qiuhe1(a,1)-qiuhe1(a,0)*qiuhe2(a,1)/(10*qiuhe2(a,0)-qiuhe1(a,0)*qiuhe1(a,0);y=a1*x+b;return y;int main()double a102=0,68,10,67.1,20,66.4,30,65.6,40,64.6,50,61.8, 60,61.0,70,60.8,80,60.4,90,60; double x,x1,q=1; char c12; int i; long n; int arw6=515,235,520,240,515,245; int arw16=315,45,320,40,325,45; int gdriver=IBM8514; int gmode=IBM8514HI; initgraph(&gdriver, &gmode, "c:TC20BGI");cleardevice();printf("input x:n");scanf("%lf",&x);printf("%fn",nihe(a,x);n=nihe(a,x)*+1;c0='y'c1='=' c4='.' for(i=10;i>1;i-) if(i!=4) ci=n%10+48; n=n/10; c11='0'x1=x;setbkcolor(7); setcolor(14);setlinestyle(0,0,3);drawpoly(3,arw);drawpoly(3,arw1);line(120,240,520,240);line(320,40,320,440);x=0;setcolor(2);setlinestyle(0,0,1); line(0+320),(int)(240-nihe(a,0)*q),(90+320),(int)(240-nihe(a,90)*q); setcolor(3); outtextxy(320,30,"Y"); outtextxy(310,245,"O"); outtextxy(525,240,"X"); outtextxy(x1+330),(240-nihe(a,x1)-10),c); settextstyle(4,0,4); outtextxy(450,400,"Nihe.");for(i=0;i<=9;i+) putpixel(ai0+320),(240-ai1*q),11);setcolor(4);setlinestyle(1,0,1);line(x1+320),(240-nihe(a,x1)*q)-80),(x1+320),(240-nihe(a,x1)*q)+120);getch();closegraph();实验结果（测试用例、实验结果）实验总结与心得通过本次实验，对曲线拟合的最小二乘法有了更深刻的了解！概念最小二乘法多项式曲线拟合，根据给定的m个点,并不要求这条曲线精确地经过这些点，而是曲线y=f(x)的近似曲线y= (x)。原理原理部分由个人根据互联网上的资料进行总结，希望对大家能有用     给定数据点pi(xi,yi)，其中i=1,2,m。求近似曲线y= (x)。并且使得近似曲线与y=f(x)的偏差最小。近似曲线在点pi处的偏差i= (xi)-y，i=1,2,.,m。 常见的曲线拟合方法:     1.使偏差绝对值之和最小          2.使偏差绝对值最大的最小          3.使偏差平方和最小          按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。推导过程：     1. 设拟合多项式为：               2. 各点到这条曲线的距离之和，即偏差平方和如下：               3. 为了求得符合条件的a值，对等式右边求ai偏导数，因而我们得到了：                                              .               4. 将等式左边进行一下化简，然后应该可以得到下面的等式：                                         .               5. 把这些等式表示成矩阵的形式，就可以得到下面的矩阵：               6. 将这个范德蒙得矩阵化简后可得到:               7. 也就是说X*A=Y，那么A = (X'*X)-1*X'*Y，便得到了系数矩阵A，同时，我们也就得到了拟合曲线。实现运行前提:1. Python运行环境与编辑环境；2. Matplotlib.pyplot图形库，可用于快速绘制2D图表，与matlab中的plot命令类似，而且用法也基本相同。代码:python  1. # coding=utf-8  2.   3. ''''' 4. 作者:Jairus Chan 5. 程序:多项式曲线拟合算法 6. '''  7. import matplotlib.pyplot as plt  8. import math  9. import numpy  10. import random  11.   12. fig = plt.figure()  13. ax = fig.add_subplot(111)  14.   15. #阶数为9阶  16. order=9  17.   18. #生成曲线上的各个点  19. x = numpy.arange(-1,1,0.02)  20. y = (a*a-1)*(a*a-1)*(a*a-1)+0.5)*numpy.sin(a*2) for a in x  21. #ax.plot(x,y,color='r',linestyle='-',marker='')  22. #,label="(a*a-1)*(a*a-1)*(a*a-1)+0.5"  23.   24. #生成的曲线上的各个点偏移一下，并放入到xa,ya中去  25. i=0  26. xa=  27. ya=  28. for xx in x:  29.     yy=yi  30.     d=float(random.randint(60,140)/100  31.     #ax.plot(xx*d,yy*d,color='m',linestyle='',marker='.')  32.     i+=1  33.     xa.append(xx*d)  34.     ya.append(yy*d)  35.   36. '''''for i in range(0,5): 37.     xx=float(random.randint(-100,100)/100 38.     yy=float(random.randint(-60,60)/100 39.     xa.append(xx) 40.     ya.append(yy)'''  41.   42. ax.plot(xa,ya,color='m',linestyle='',marker='.')  43.   44.   45. #进行曲线拟合  46. matA=  47. for i in range(0,order+1):  48.     matA1=  49.     for j in range(0,order+1):  50.         tx=0.0  51.         for k in range(0,len(xa):  52.             dx=1.0  53.             for l in range(0,j+i):  54.                 dx=dx*xak  55.             tx+=dx  56.         matA1.append(tx)  57.     matA.append(matA1)  58.   59. #print(len(xa)  60. #print(matA00)  61. matA=numpy.array(matA)  62.   63. matB=  64. for i in range(0,order+1):  65.     ty=0.0  66.     for k in range(0,len(xa):  67.         dy=1.0  68.         for l in range(0,i):  69.             dy=dy*xak  70.         ty+=yak*dy  71.     matB.append(ty)  72.    73. matB=numpy.array(matB)  74.   75. matAA=numpy.linalg.solve(matA,matB)  76.   77. #画出拟合后的曲线  78. #print(matAA)  79. xxa= numpy.arange(-1,1.06,0.01)  80. yya=  81. for i in range(0,len(xxa):  82.     yy=0.0  83.     for j in range(0,order+1):  84.         dy=1.0  85.         for k in range(0,j):  86.             dy*=xxai  87.         dy*=matAAj  88.         yy+=dy  89.     yya.append(yy)  90. ax.plot(xxa,yya,color='g',linestyle='-',marker='')  91.   92. ax.legend()  93. plt.show()  专心-专注-专业