数值分析实验报告（拉格朗日插值牛顿插值最小二乘法）
2010-06-02 18:33阅读：
实验1 拉格朗日插值法
一、方法原理
n次拉格朗日插值多项式为：Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x)
n=1时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+
y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0)
n=2时，称为二次插值或抛物线插值，精度相对高些
L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1)
二、主要思路
使用线性方程组求系数构造插值公式相对复杂，可改用构造方法来插值。
对节点xi(i=0,1,…,n)中任一点xk(0<=k<=n)作一n
次多项式lk(xk)，使它在该点上取
值为1，而在其余点x
i(i=0,1,…,k-1,k+1,…,n)上为0，则插值多项式为L
n(x)=y
0l
0(x)+y
1l
1(x)+y
2l
2(x)+…+y
nl
n(x)
上式表明：n
个点x
i(i=0,1,…,k-1,k+1,…,n)都是l
k(x)的零点。可求得l
k
三．计算方法及过程：1.输入节点的个数n
2.输入各个节点的横纵坐标
3.输入插值点
4.调用函数，返回z
函数语句与形参说明
程序源代码如下：
形参与函数类型
|
参数意义
|
int n
|
节点的个数
|
double x[n]（double *x）
|
存放n个节点的值
|
double y[n]（double *y）
|
存放n个节点相对应的函数值
|
double p
|
指定插值点的值
|
double fun()
|
函数返回一个双精度实型函数值，即插值点p处的近似函数值
|
#include<iostream>
#include<math.h>
using namespace std;
#define N 100
double fun(double *x,double *y, int n,double p);
void main()
{int i,n;
cout<<'输入节点的个数n：';
cin>>n;
double x[N], y[N],p;
cout<<'please input xiangliang x= '<<endl;
for(i=0;i<n;i++)cin>>x[i];
cout<<'please input xiangliang y= '<<endl;
for(i=0;i<n;i++)cin>>y[i];
cout<<'please input LagelangrichazhiJieDian p=
'<<endl;
cin>>p;
cout<<'The Answer= '<<fun(x,y,n,p)<<endl;
system('pause') ;}
double fun(double x[],double y[], int n,double p)
{double z=0,s=1.0;
int k=0,i=0;
double L[N];
while(k<n)
{ if(k==0)
{ for(i=1;i<n;i++)s=s*(p-x[i])/(x[0]-x[i]);
L[0]=s*y[0];
k=k+1;}
else
{s=1.0;
for(i=0;i<=k-1;i++)s=s*((p-x[i])/(x[k]-x[i]));
for(i=k+1;i<n;i++)
s=s*((p-x[i])/(x[k]-x[i]));
L[k]=s*y[k];
k++;}
}
for(i=0;i<n;i++)z=z+L[i];
return z;
}
四．运行结果测试：
五．实验分析
n=2时，为一次插值，即线性插值
n=3时，为二次插值，即抛物线插值
n=1，此时只有一个节点，插值点的值就是该节点的函数值
n<1时，结果都是返回0的；这里做了n=0和n=-7两种情况
3<n<100时，也都有相应的答案
常用的是线性插值和抛物线插值，显然，抛物线精度相对高些
n次插值多项式Ln（x）通常是次数为n的多项式，特殊情况可能次数小于n.例如：通过三点的二次插值多项式L2(x),如果三点共线，则y=L2(x)就是一条直线，而不是抛物线，这时L2(x)是一次式。
拟合曲线光顺性差
实验2 牛顿插值法
一、方法原理及基本思路
在拉格朗日插值方法中，若增加一个节点数据，其插值的多项式需重新计算。现构造一个插值多项式Nn(x)，只需对Nn-1(x)作简单修正(如增加某项)即可得到，这样计算方便。
利用牛顿插值公式，当增加一个节点时，只需在后面多计算一项，而前面的计算仍有用；另一方面Nn(x)的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。
由线性代数知，对任何一个不高n次的多项式P(x)=b0＋b1x＋b2x2＋…＋bnxn (幂基) ①
也可将其写成P(x)=a
0＋a
1(x－x
0)＋a
2(x－x
0)
(x－x
1)＋…＋a
n(x－x
0)
…(x－x
n-1)
其中ai为系数，xi为给定节点，可由①求出ai
一般情况下，牛顿插值多项式Nn(x)可写成：
Nn(x)=
a
0＋a
1(x－x
0)＋a
2(x－x
0)
(x－x
1)＋…＋a
n(x－x
0)
…(x－x
n-1))
只需求出系数ai，即可得到插值多项式。
二、计算方法及过程
1.先后输入节点个数n和节点的横纵坐标，插值点的横坐标，最后输入精度e
2. 用do-while循环语句得到跳出循环时k的值
3.将k值与n-1进行比较，若在达到精度时k<n-1,则输出计算结果；若此时k=n-1,则计算失败！
函数语句与形参说明
形参与函数类型
|
参数意义
|
int n
|
节点的个数
|
float x[MAX]
|
存放n个节点的值（从小到大）
|
Float y[MAX];
|
存放n个节点相对应的函数值
|
float x0,y0
|
指定插值点的横纵坐标
|
float e
|
精度
|
程序源代码如下：
#include<iostream>
#include<math.h>
using namespace std;
#define MAX 100
void main()
{
float x[MAX],y[MAX];
float x0,y0,e,N1;
float N0=0;
int i,k,n;
cout<<'输入节点的个数,n=';
cin>>n;
cout<<'请输入插值节点！'<<endl;
for(i=0;i<n;i++) cin>>x[i];
cout<<'请输入对应的函数值！'<<endl;
for(i=0;i<n;i++) cin>>y[i];
cout<<'请输入插值点,x0=';
cin>>x0;
cout<<'请输入精度,e=';
cin>>e;
y0=1;
N1=y[0];
k=0;
do
{
k++;
N0=N1;
y0=y0*(x0-x[k-1]);
for(i=0;i<k;i++) y[k]=(y[k]-y[i])/(x[k]-x[i]);
N1=N0+y0*y[k];
}
while (fabs(N1-N0) > e && k<(n-1));
if (k==(n-1)) cout<<'计算失败！';
if (k<(n-1))
cout<<'输出结果y[x0]='<<N1<<endl;
system('pause');
}
三．运行结果测试：
题一：已知f(x)=sh(x)的函数表如下：计算f(0.23)的近似值
xi
|
0
|
0.20
|
0.30
|
0.50
|
Sh(xi)
|
0
|
0.20134
|
0.30452
|
0.52110
|
题二：已知根号100等于10，根号121等于11，根号144等于12；
用二次插值和三次插值计算根号115的近似值。
四．实验分析
显然，题二中当n=2,n=3时，k的值没有在n-1前满足各自的精度要求，所以计算失败
从题一中可以看出，插值点的个数、精度、插值点的选择都会影响实验的结果；
我们通常会选择与插值点最接近的节点，可以提高精度；
在可以计算出结果的情况下，插值点越多，结果越精确。
实验3 最小二乘法求拟合曲线
一、方法原理及思路
已知数据对(xj, yj)(j=1,2,…,n)，求多项式
为最小，这就是一个最小二乘问题。
二、计算方法及过程：
1.
先后输入数据的对数m,近似多项式的最高次数n及各个数据的横纵坐标
2.
用for循环语句得到线性方程组Ga=d
3.
用高斯列主元消去法得到的解即为近似多项式的系数ai(i=0,1,…,n)
函数语句与形参说明
形参与函数类型
|
参数意义
|
int m
|
数据的对数,即坐标的个数
|
int n
|
近似多项式的最高次数
|
float x[MAX]
|
存放m个横坐标的值（从小到大）
|
Float y[MAX];
|
存放m个纵坐标的值（与相应的横坐标对应）
|
Float a[ MAX +1],
|
近似多项式的系数
|
float G[ MAX +1][ MAX +1]
|
法方程的系数矩阵，且为对称的，正定的
|
float s[2 * MAX +1]
|
系数矩阵G中的元素
|
Float b[ MAX +1];
|
法方程中等号右边的列向量
|
程序原代码如下
#include<math.h>
#include<iostream>
#define
MAX
100
using namespace std;
int main()
{ float t;
int i,j,k,r,m,n;
cout<<'请输入数据的对数，m=';
cin>>m;
cout<<'请输入近似多项式的最高次数，n=';
cin>>n;
float y[ MAX],x[ MAX ],x1[ MAX ],a[ MAX +1],G[ MAX
+1][ MAX +1], s[2* MAX +1],b[ MAX +1];
cout<<'请输入各个数据的横坐标！'<<endl;
for (i=0;i<m;i++)cin>>x[i];
cout<<'请输入各个数据的纵坐标！'<<endl;
for
(i=0;i<m;i++)cin>>y[i];
s[0]=m; b[0]=0.0;
for(j=0;j<m;j++)
{b[0]=b[0]+y[j];
x1[j]=1.0;
}
for (i=1;i<=n;i++)
{b[i]=0.0; s[2*i-1]=0.0;s[2*i]=0.0;
for(j=0;j<m;j++)
{x1[j]=x1[j]*x[j];
b[i]=b[i]+x1[j]*y[j];
s[2*i-1]=s[2*i-1]+ x1[j]*x1[j]/x[j];
s[2*i]=s[2*i]+x1[j]*x1[j];}
}
for(i=0;i<=n;i++)
{ for(j=0;j<=n;j++) G[i][j]=s[i+j]; }
for(k=0;k<n;k++)
//高斯列主元消去法
{
t=abs(G[k][k]);
r=k;
for(i=k+1;i<=n;i++)
{
if(t<abs(G[i][k]))//选主元
{
t=G[i][k];
r=i; }
}
if(r>k)//行交换
{
for(j=k;j<=n;j++)
{
t=G[k][j];
G[k][j]=G[r][j];
G[r][j]=t; }
t=b[k];
b[k]=b[r];
b[r]=t; }
for(i=k+1;i<=n;i++)//消元
{
t=G[i][k]/G[k][k];
for(j=k+1;j<=n;j++) G[i][j]=G[i][j]-t*G[k][j];
b[i]=b[i]-t*b[k]; }
}
a[n]=b[n]/G[n][n];//回代
for(k=1;k<=n;k++)
{
i=n-k; a[i]=b[i];
for(j=i+1;j<=n;j++)a[i]=a[i]-G[i][j]*a[j];
a[i]=a[i]/G[i][i]; }
cout<<'输出结果a'<<endl;
//输出答案
for(i=0;i<=n;i++)cout<<a[i]<<'\t';
system('pause');
}
三：运行结果测试：
四．实验分析
由于曲线拟合的最小二乘法一般是经过描点，确定其近似多项式的形式，但由于给定的点有误差，所以拟合曲线的数学模型并不是一开始就能选的好的，往往要通过分析确定若干模型后，再经过实际计算，才能选到较好的模型。
对我的博文感兴趣的话，关注我哦~~~~
