最小二乘法和Romberg求积分.doc

计算方法实验报告

湖北师范学院计算机科学与技术学院

一、实验题目
最小二乘法和Romberg求积分
二、实验内容
1、最小二乘法拟合
本实验是对于已知的m+1对不带权系数的离散数据，计。在连续的函数空间中选定n+1个线性无关的基函数来构造拟合函数，求得当系数是最优化模型时的解。
2、romberg求积分法
Romberg算法是利用复合梯形公式，在对积分区间的步长逐次折半的过程中，求得积分I=∫abf(x)dx的近似值序列{T2k}。

三、法思想描述
1、 最小二乘法拟合
（1）、根据最小二乘法的拟合原理将问题简化为求由待拟合离散数据的正则方程组的求解问题；
（2）、输入一散数据，存放于数组中，并给定拟合次数n根据待拟合的次数确定对应的系数矩阵为范德蒙行列式的线性方程组；
（3）、由得到正则方程组的系数矩阵，如下式，同时由得到值向量；

Y1
Y2

Ym

＝
a0
a1
a2
…
an
1 x0 x0 2… x0n
1 x1 xin… xin
… 　…　　…　 …
…　　…　　…　…
1 xm xm2… xmn

(4)、用最列主元素法解出正则方程组的根，就得到“二”中描述的系数的最优解。
2、romberg求积分法
Romberg算法是区间主次二分的过程中，对梯形值进行加权平均以获得准确程度较高的积分值的一种方法。
对于定积分I= 的复化梯形公式 ，其余项
将积分区间[a，b]，逐次折半，假设 ,以保证复化梯形公式余项系数是非零的，则构成相应的外推法称为Romberg算法：
下标m外推得到的第m个算法。 中的求和项包括了每次外推后新增加点上的函数值。注意对某个k，被积函数有性质 。说明余项中 ，则(*)式中的外推法要作相应的修改，否则外推可能失效。

四、设计说明

五、源程序清单
#include<iostream.h>
#include<math.h>
#include<stdio.h>
#define N 15
#define eps 0.00000001
int n,k;
double X[50],Y[50],A[50][50],F[50][50],B[50],root[50];
double F1[N],F2[N];
int m;
int v=0;
double aa,bb;
void getmessage(){
cout<<"输入待拟合数据的组数:";
cin>>n;
cout<<"输入待拟合数据的次数：";
cin>>k;
while(n<k){
cout<<"拟合次数输入过大，请重新输入：";
cin>>k;
}
cout<<"输入待拟合数据x,y:"<<endl;
for(int i=0;i<n;i++){
cin>>X[i];
cin>>Y[i];
}
cout<<"待拟合数据输入完毕！"<<endl;
}
void getfar(){
int i,j,t;
double s=0;
for(i=0;i<n;i++)
for(j=0;j<=k;j++)
A[i][j]=pow(X[i],j);
for( t=0;t<=k;t++)
for(i=0;i<=k;i++){
for(j=0;j<n;j++)
s=s+A[j][t]*A[j][i];
F[t][i]=s;
s=0;
}
for(i=0;i<=k;i++){
for(j=0;j<n;j++)
s=s+A[j][i]*Y[j];
B[i]=s;
s=0;
}
for(i=0;i<=k;i++){
for(j=0;j<=k;j++)
cout<<F[i][j]<<" ";
cout<<endl;
}
for(j=0;j<=k;j++)
cout<<"b"<<j<<"="<<B[j]<<endl;
}

void liezhu(){
double tt,s=0,f1;
int i,j,t,max;
for(i=0;i<=k;i++){
max=i;
for(j=i+1;j<=k;j++)
if(F[j][i]>F[max][i]) max=j;
if(max!=i){
for(t=i;t<=k;t++){
tt=F[i][t];
F[i][t]=F[max][t];
F[max][t]=tt;
}
tt=B[i];
B[i]=B[max];
B[max]=tt;
}
for(j=i+1;j<=k;j++){
f1=F[j][i]/F[i][i];
for(t=i;t<=k;t++)
F[j][t]=F[j][t]-f1*F[i][t];
B[j]=B[j]-f1*B[i];
}
}
root[k]=B[k]/F[k][k];
for(i=k-1;i>=0;i--){
for(j=k;j>i;j--)
s=s+F[i][j]*root[j];
root[i]=(B[i]-s)/F[i][i];
s=0;
}
for(i=0;i<=k;i++){
cout<<"a"<<i<<"="<<root[i]<<",";
}
cout<<endl;
}

void ercheng(){
getmessage();
getfar();
liezhu();
}

void Romberg()
{
int k=0,n=1,d,i,j;
float a,b,h;
double f(float x);
double y,c,T[50];
printf("**********Romberg求积算法**********\n\n\n");
printf("请输入下限a:a=");
scanf("%f",&a);
printf("\n");
printf("请输入上限b:b=");
scanf("%f",&b);
printf("\n");
printf("\n\n\n");
h=b-a;
T[0]=h/2*(f(a)+f(b));
printf("计算过程T-数表如下：\n");
printf("%.8lf\n",T[0]);
c=T[0];
aaa: k=k+1;
y=0;
for(i=1;i<=n;i++)
{
y=y+f(a+(i-0.5)*h);
}
T[k]=(T[k-1]+h*y)/2;
printf("%.8lf ",T[k]);
d=1;
for(j=k;j>=1;j--)
{
d=4*d;
T[j-1]=T[j]+(T[j]-T[j-1])/(d-1);
printf("%.8lf ",T[j-1]);
}
printf("\n");
if(fabs(c-T[0])<eps)
{
printf("\n\t\t满足精度要求的值I≈%.8lf\n\n",T[0]);
printf("\t\t误差值为ε=%.8le\n",fabs(c-T[0]));
}
else
{
n=2*n;
h=h/2;
c=T[0];
goto aaa;
}
}
double f(float x)
{
double z;
if(x==0)
z=1;
else
z=sin(x)/x;
//z=sqrt(x);
//z=4/(1+pow(x,2));
return z;
}

void main(){
int ch;
char yn;
while(1)
{
cout<<"/**** 功能选择 ****/"<<endl;
cout<<"/**** 1.最小二乘拟合 ****/"<<endl;
cout<<"/**** 2.Romberg算法 ****/"<<endl;
cout<<"请输入:";
cin>>ch;
switch(ch)
{
case 1:ercheng();
break;
case 2:Romberg();
break;
default:
cout<<"退出系统吗?(y/n):";
cin>>yn;
if(yn=='Y'||yn=='y'){
cout<<"谢谢使用!"<<endl;
return;
}
else
break;
}
}
cout<<"是否继续(y/n):";
cin>>yn;
if(yn!='Y'&&yn!='y')
return;
}
六、测试结果
最小二乘拟合测试结果如下：

Romberg 测试结果：