《计算方法》实验报告学院:计算机学院专业:计算机科学与技术指导教师:JW-++1 爨莹班级学号:201207010229 姓名:图尔荪托合提实验三曲线拟合的最小二乘法1、实验目的:在科学研究与工程技术中,常常需要从一组测量数据出发,寻找变量的函数关系的近似表达式,使得逼近函数从总体上与已知函数的偏差按某种方法度量能达到最小而又不一定过全部的点。
这是工程中引入最小二曲线拟合法的出发点。
充分掌握:1.最小二乘法的基本原理;2.用多项式作最小二乘曲线拟合原理的基础上,通过编程实现一组实验数据的最小二乘拟合曲线。
2、实验要求:1) 认真分析题目的条件和要求,复习相关的理论知识,选择适当的解决方案和算法;2) 编写上机实验程序,作好上机前的准备工作;3) 上机调试程序,并试算各种方案,记录计算的结果(包括必要的中间结果);4) 分析和解释计算结果;5) 按照要求书写实验报告;3、实验内容:1) 给定数据如下:x :0.15,0.4,0.6 ,1.01 ,1.5 ,2.2 ,2.4,2.7,2.9,3.5 ,3.8 ,4.4,4.6 ,5.1 ,6.6,7.6;y :4.4964,5.1284,5.6931 ,6.2884 ,7.0989 ,7.5507 ,7.5106,8.0756,7.8708,8.2403 ,8.5303 ,8.7394,8.9981 ,9.1450 ,9.5070,9.9115;试作出幂函数拟合数据。
2) 已知一组数据:x :0,0.1,0.2 ,0.3 ,0.4 ,0.5 ,0.6,0.7,0.8,0.9 ,1y :-0.447,1.978,3.28 ,6.16 ,7.08 ,7.34 ,7.66,9.56,9.48,9.30 ,11.2;试用最小二乘法求多项式函数,使与此组数据相拟合。
4、题目:曲线拟合的最小二乘法5、原理:从整体上考虑近似函数同所给数据点(i=o,1,…,m误差(i=o,1,…,m)的大小,常用的方法有以下三种:一是误差(i=0,1,…口绝对值的最大值,即误差向量的%—范数;二是误差绝对值的和,即误差向量r的1 —范数;三是误差平方和的算术平方根,即误差向量r的2—范数;前两种方法简单、自然,但不便于微分运算,后一种方法相当于考虑2—范数的平方,因此在曲线拟常采用误差平方和来度量误差(i=0,1,…,m)的整体大小.。
数据拟合的具体作法是:对给定数据(i=0,1,…m),在取定的函数类中,求,使误差(i=0,1,…,m)的平方和最小。
6设计思想:从几何意义上讲,就是寻求与给定点(i=0,1,…,m)的距离平方和为最小的曲线。
函数称为拟合函数或最小二乘解,求拟合函数的方法称为曲线拟合的最小二乘法。
7、对应程序:(1)幕函数程序#i nclude<stdio.h>#in clude<math.h>void mai n(){ double a0[16][2],a1[2][16],A[2][2],Y[2];int i,j,k;double x[16]={0.15,0.4,0.6,1.01,1.5,2.2,2.4,2.7,2.9,3.5,3.8,4.4,4.6,5.1,6.6, 7.6},y[16]={4.4964,5.1284,5.6931,6.2884,7.0989,7.5507,7.5106,8.0756,7.8708,8.2403,8.5303,8.7394,8.9981,9.1450,9.5070,9.9115};double m0,m1, n;for(i=0;i<=15;i++){a0[i][0]=1;a0[i][1]=log(x[i]);y[i]=log(y[i]);printf(" 输入X 的值:");prin tf("%f ",a0[i][1]);prin tf(" 得到的对应的函数值:");prin tf("%f \n ",y[i]);}prin tf("\n");for(i=0;i<=15;i++) for(j=0;j<=1;j++){a1[j][i]=a0[i][j];}// 以上正确A[0][0]=0;A[0][1]=0;A[1][0]=0;A[1][1]=0;Y[0]=0;Y[1]=0;for(i=0;i<=1;i++){for(j=0;j<=1;j++){for(k=0;k<=15;k++){A[i][j]+=a1[i][k]*a0[k][j];}}for(i=0;i<=1;i++) {for(j=0;j<=15;j++){Y[i]+=a1[i][j]*y[j];}} m0=(Y[0]*A[1][1]-Y[1]*A[0][1])/(A[0][0]*A[1][1]-A[1][0]*A[0][1]); n=(Y[0]*A[1][0]-Y[1]*A[0][0])/(A[0][1]*A[1][0]-A[1][1]*A[0][0]);m1=exp(m0);printf(" 得到的幂函数X 的系数是:%f\n",m1); printf(" 得到的幂函数X 的指数是: %f\n",n);} (2)最小二乘法求多项式#include <stdio.h> #include <conio.h> #include <math.h> #include <process.h> #define N 11//N 个点#define T 3 //T 次拟合#define W 1// 权函数#define PRECISION 0.00001 float pow_n(float a,int n) {int i;if(n==0) return(1);float res=a; for(i=1;i<n;i++) {res*=a;} return(res);} void mutiple(float a[][N],floatb[][T+1],float c[][T+1]) {float res=0;int i,j,k; for(i=0;i<T+1;i++) for(j=0;j<T+1;j++){res=0; for(k=0;k<N;k++) {res+=a[i][k]*b[k][j];c[i][j]=res;}}}void matrix_trans(float a[][T+1],float b[][N]){int i,j; for(i=0;i<N;i++) {for(j=0;j<T+1;j++) {b[j][i]=a[i][j];}}}void init(float x_y[][2],int n){int i;prin tf(" 请输入%c个已知点:\n",N);for(i=0;i<n;i++){printf("(x%d y%d):",i,i);scanf("%f %f",&x_y[i][0],&x_y[i][1]);}}void get_A(float matrix_A[][T+1],float x_y[][2],int n){int i,j;for(i=0;i<N;i++){for(j=0;j<T+1;j++) {matrix_A[i][j]=W*pow_n(x_y[i][0],j);}} }void print_array(float array[][T+1],int n){int i,j; for(i=0;i<n;i++) {for(j=0;j<T+1;j++) {printf("%-g ",array[i][j]);} printf("\n");}}void convert(float argu[][T+2],int n){int i,j,k,p,t;float rate,temp;for(i=1;i<n;i++){for(j=i;j<n;j++) {if(argu[i-1][i-1]==0) {for(p=i;p<n;p++){if(argu[p][i-1]!=0)break;}if(p==n){printf(" 方程组无解!\n");exit(0);} for(t=0;t<n+1;t++){temp=argu[i-1][t];argu[i-1][t]=argu[p][t]; argu[p][t]=temp;} }rate=argu[j][i-1]/argu[i-1][i-1]; for(k=i-1;k<n+1;k++) {argu[j][k]-=argu[i-1][k]*rate;if(fabs(argu[j][k])<=PRECISION) argu[j][k]=0;}}}}void compute(float argu[][T+2],int n,float root[]){int i,j;float temp;for(i=n-1;i>=0;i--){temp=argu[i][n];for(j=n-1;j>i;j--){temp-=argu[i][j]*root[j];}root[i]=temp/argu[i][i];}}void get_y(float trans_A[][N],float x_y[][2],float y[],int n){int i,j;float temp;for(i=0;i<n;i++){temp=0;for(j=0;j<N;j++){temp+=trans_A[i][j]*x_y[j][1];}y[i]=temp;}}void cons_formula(float coef_A[][T+1],float y[],float coef_form[][T+2]) {int i,j; for(i=0;i<T+1;i++){for(j=0;j<T+2;j++){if(j==T+1)coef_form[i][j]=y[i];elsecoef_form[i][j]=coef_A[i][j];}}}void print_root(float a[],int n){int i,j;printf("%d 个点的%4次拟合的多项式系数为:\n",N,T);for(i=0;i<n;i++){printf("a[%d]=%g,",i+1,a[i]);}printf("\n");printf(" 拟合曲线方程为:\ny(x)=%g",a[0]);for(i=1;i<n;i++){printf(" + %g",a[i]);for(j=0;j<i;j++){printf("*X");}}printf("\n");}void process(){float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_form u[T+1][T+2],y[T+1],a[T+1];init(x_y,N);get_A(matrix_A,x_y,N);"个点的2次拟合的多顶式系数为=丑【11=-0・510421「日[2]=29・3249,a[3 1=-33 *5252 .a [4]=16 .0762, 执合曲线方程为:</<x>=-0,610421 ■* 29.3249**« - -33.9252*X*K * 16,6762«X*X*K匚也Lk e r^\Admiriirtrater\Des Ictc p\D 启bug'吉吕m a.exe-"-1.897120 -9.91t29i -0.510S2C 0 .価 950 0.405465 0.788457 缜8?54&? 8.993252 1.064711 1.252763 1.3350B1 1.481605 1.52&056 1.&29241 1.887S70 2,020148L.S03277 i.&34794 1.739255 1.838707 L ・?5?940 Z.021640 2.016315 Z ・6SSS47&3160 2.109037 Z.143&25 2.1&7S42 2.197013 2.21320P 2.2&2028 2.293696jpi| j^nl j^Irl j^TTI数数应应应应需T 1—T叫T r亠、■•--,-『L~^TQ 「T二rrEl-D一・「1■『「■E_i “Jboi寸曰寸目6.4157460.209345'rcss anu ]«ey to continu.e —^n-A,Hr __i^n LLJlf*rf —rtnJ^o一-■气『■J^n-u A y A y A y A y A y A y A .A ^Ay A^(2)最小二乘法求多项式结果10、实验体会:最小二乘法虽然看起来简单,但它在数值计算及应用上却非常重要,通过C语言程序编写及运行最小二乘法的程序,加强了我们的综合能力,这些告诉我们实践很重要,实验中我们最应该注意的还是细心,认真。