软件下载 | 资讯教程 | 最近更新 | 下载排行 | 一键转帖 | 发布投稿
您的位置:最火下载站 > 电脑教程 > 编程开发 > C/C++开发 > 数值计算:插值的数值解法的C++程序

数值计算:插值的数值解法的C++程序

插值的各种算法的C++实现,是学数值计算方法的实验,有拉格朗日插值,牛顿插值和调用gsl函数进行的样条插值。 程序使用了gsl以及我自己写的两个库(均已上传)。

Linux下编译方法:g++ interp.cpp -o interp -lm -lgsl -lgslcblas

附图为把程序对y=sin(x)在[1,20]区间内插值运行结果复制到scilab里画的图,实线是原始函数图像,*为拉格朗日插值结果,虚线为牛顿插值结果,点实线为gsl函数插值的结果。

代码:

//数值计算实验 插值

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <mylib/MyClock.h>
#include <mylib/InitData.h>

using namespace std;

//要插值的函数,为sin(x),可改。
double f(double x)
{
return sin(x);
}

//产生函数区间为[a,b]上的n个点的插值函数数值,x,y分别存入横纵座标
bool InitFun(const double & a, const double & b, const int & n,
vector<double> & x, vector<double> & y)
{
//如果a大于b返回错误。
if (a > b)
{
return false;
}
//计算x值间隔大小,将[a,b]均分为n-1份,这样连端点,一共n个点
double gap = (b-a)/(n-1.0);
x.resize(n);
y.resize(n);
for (int i = 0; i < n; i++)
{
x[i] = a + gap*i;
y[i] = f(x[i]);
}
return true;
}

//根据x,y,n,计算n阶拉格朗日插值y1 = f(x1)
bool Lagrange(const vector<double> & x, const vector<double> & y,
const double & x1, double & y1)
{
y1 = 0.0;
int n = x.size();
for (int j = 0; j < n; j++)
{
double L = 1.0;
for (int i = 0; i < n; i++)
{
if (i != j)
{
L *= (x1-x[i])/(x[j] - x[i]);
}
}
y1 += L*y[j];
}
return true;
}

//根据x,y,n,计算n阶牛顿插值y1 = f(x1),
//参考http://hi.baidu.com/trady830/blog/item/6745631e4ab4366af624e4a8.html
bool Newton(const vector<double> & x, const vector<double> & y,
const double & x1, double & y1)
{
int n = x.size();

//计算差商
vector<double> m_diff(n);

//先存放0阶差商
for (int i = 0; i < n; i++)
{
m_diff[i] = y[i];
}
for (int i = 0; i < n-1; i++)
{
for (int j = n-1; j > i; j--)
{
m_diff[j] = (m_diff[j] - m_diff[j-1])/(x[j] - x[j-1-i]);
}
}
double temp = 1.0;
y1 = m_diff[0];
for (int i = 0; i < n-1; i++)
{
temp = temp * (x1 - x[i]);
y1 += temp*m_diff[i+1];
}
return true;
}

//根据x,y,n,调用gsl函数计算n阶插值y1 = f(x1),
//参考gsl文档
bool GSLInterp(vector<double> & x, vector<double> & y,
double & x1, double & y1)
{
int n = x.size();
double *gx, *gy;
try
{
gx = new double[n];
gy = new double[n];
}
catch (bad_alloc)
{
cout<<"内存分配失败!"<<endl;
return false;
}

for (int i = 0; i < n; i++)
{
gx[i] = x[i];
gy[i] = y[i];
}

gsl_interp_accel *acc = gsl_interp_accel_alloc();
gsl_spline *spline = gsl_spline_alloc(gsl_interp_cspline, n);
gsl_spline_init(spline, gx, gy, n);

y1 = gsl_spline_eval(spline, x1, acc);

gsl_spline_free(spline);
gsl_interp_accel_free(acc);
delete[] gx;
delete[] gy;
return true;
}

int main()
{
int n;
double a, b;
cout<<"请输入插值区间:";
cin>>a>>b;
cout<<"请输入样本量:";
cin>>n;
if (a > b || n <= 0)
{
cout<<"输入错误!"<<endl;
exit(1);
}
vector<double> x(n);
vector<double> y(n);

//产生插值数据
if (!InitFun(a, b, n, x, y))
{
cout<<"产生插值数据出错!"<<endl;
exit(1);
}
/* //调试InitFun函数用的
else
{
cout<<"产生的插值数据为:"<<endl;
cout.precision(10);
cout.width(12);
cout<<"x = "<<endl;
ShowData(x, cout);
cout<<"y = "<<endl;
ShowData(y, cout);
}
*/
//设定输出精度
cout.precision(8);
cout.width(10);

//测试插值函数
/*
double x1, y1;
cout<<"输入插值区间["<<a<<","<<b<<"]内要插值的x值:";
cin>>x1;
if (x1 < a || x1 > b)
{
cout<<"输入错误!"<<endl;
exit(1);
}
//if (Lagrange(x, y, x1, y1))
//if (Newton(x, y, x1, y1))
if (GSLInterp(x, y, x1, y1))
{
//cout<<"拉格朗日插值结果为:"<<y1<<endl;
//cout<<"牛顿插值结果为:"<<y1<<endl;
cout<<"gsl函数插值结果为:"<<y1<<endl;
//cout<<"所求的精确值为:"<<sqrt(abs(x1))<<endl;
cout<<"所求的精确值为:"<<sin(x1)<<endl;
}
else
{
cout<<"插值失败!"<<endl;
}
*/

//批量插值
int m;
cout<<"输入要插值的数据数量:";
cin>>m;
if (m <= 0)
{
cout<<"输入数值必须大于0"<<endl;
exit(1);
}
vector<double> x2(m);
vector<double> y2(m); //拉格朗日插值结果
vector<double> y3(m); //牛顿插值结果
vector<double> y4(m); //gsl插值结果
double gap = (b-a)/(m-1);
for (int i = 0; i < m; i++)
{
x2[i] = a + gap*i;
if (Lagrange(x, y, x2[i], y2[i]) != true)
{
cout<<"拉格朗日插值失败,程序退出!"<<endl;
exit(1);
}
if (Newton(x, y, x2[i], y3[i]) != true)
{
cout<<"牛顿插值失败,程序退出!"<<endl;
exit(1);
}
if (GSLInterp(x, y, x2[i], y4[i]) != true)
{
cout<<"gsl插值失败,程序退出!"<<endl;
exit(1);
}
}
cout<<"插值情况为:"<<endl;
cout<<"x = "<<endl;
ShowData(x2, cout);
cout<<"拉格朗日插值结果:"<<endl;
cout<<"y = "<<endl;
ShowData(y2, cout);
cout<<"牛顿插值结果:"<<endl;
cout<<"y = "<<endl;
ShowData(y3, cout);
cout<<"gsl插值结果:"<<endl;
cout<<"y = "<<endl;
ShowData(y4, cout);
return 0;
}

    相关阅读
    栏目导航
    推荐软件