CPP&MATLAB实现拉格朗日插值法

Posted on 2016-09-05 10:10  不秩稚童  阅读(4122)  评论(1编辑  收藏  举报

    开始学习MATLAB(R和Python先放一放。。。),老师推荐一本书,看完基础就是各种算法。。。首先是各种插值。先说拉格朗日插值法,这原理楼主完全不懂的,查的维基百科,好久才看懂。那里讲的很详细,这里就不在赘述了。一般看这个范例,在回头看公式就比较容易理解。

 

关于MATLAB的实现,查了很多资料,下面的版本最好理解。

%lagran1.m 
%求拉格朗日插值多项式和基函数 
%输入的量:n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量X,纵坐标向量Y 
%输出的量:n次拉格朗日插值多项式L,基函数l
function [L,l] = lagran1(X,Y)
% m为已知插值点的个数
m = length(X); 
for k = 1 : m
    % 迭代相乘前的初始化
    % V代表每个k下的拉格朗日基本多项式
    V = 1; 
    for i = 1 : m 
        if k ~= i 
            % conv这里是用于多项式相乘
            % poly(a) 当a为数时,返回多项式(x - a)
            V = conv(V,poly(X(i))) / (X(k) - X(i)); 
        end 
    end
    %l将所有的拉格朗日多项式存为一列
    %poly2sym就是转化为多项式,下文有简单的demo
    l(k, :) = poly2sym(V); 
end 
% 得到最终的拉格朗日多项式
L = Y * l;

  关于poly2sym,是将数组转化为多项式的函数,不指定符号,默认为x.

   运行函数:

  查看拉格朗日多项式:

  运用此函数进行插值:

  接下来,最近学了一点cpp,就想用用cpp实现一遍,查到了这篇博文。这里就直接贴上了。

#include<iostream>
#include<string>
#include<vector>
using namespace std;
/*
vector:不是一种数据类型。 vector<int>是一种数据类型。
vector容器是一个模板类,可以存放任何类型的对象(但必须是同一类对象)。
vector对象可以在运行时高效地添加元素,并且vector中元素是连续存储的。
*/

// 关于 & 的解释
/*
当我们声明一个变量的同时,它必须被存储到内存中一个具体的单元中。
通常我们并不会指定变量被存储到哪个具体的单元中—幸亏这通常是由编译器和操作系统自动完成的,
但一旦操作系统指定了一个地址,有些时候我们可能会想知道变量被存储在哪里了。
这可以通过在变量标识前面加与符号ampersand sign (&)来实现,它表示"...的地址" ("address of"),
因此称为地址操作符(adress operator),又称去引操作符(dereference operator)
*/

// 函数声明,是的其在完整定义前可以被(main函数)使用
double Lagrange(int N,vector<double>&X,vector<double>&Y,double x);

int main()
{
	char a='n';
	do{
		cout<<"请输入插值次数n的值:"<<endl;
		int N;
		cin>>N;
		vector<double>X(N,0);
		vector<double>Y(N,0);
		cout<<"请输入插值点对应的值及其函数值(Xi,Yi):"<<endl;
		for(int a=0;a<N;a++){
			cin>>X[a]>>Y[a];
		}
		cout<<"请输入要求值的x的值:"<<endl;
		double x;
		cin>>x;
		double result=Lagrange(N,X,Y,x);
		cout<<"由拉格朗日插值法得出的结果:"<<result<<endl;
		cout<<"是否要继续?<y/n>:";
		cin>>a;
	}while(a=='y');
	return 0;
}

//vector<double>,创建空容器,其对象类型为double类
double Lagrange(int N,vector<double>&X,vector<double>&Y,double x)
{
	double result=0;
	for(int i=0;i<N;i++){
		double temp=Y[i];
		for(int j=0;j<N;j++){
			if(i != j){
				temp = temp*(x-X[j]);
				temp = temp/(X[i]-X[j]);
			}
		}// 与上文MATLAB实现不同,这里直接将拉格朗日多项式的各项直接相加
		result += temp;
	}
	return result; 
}

  关于 & 的用法:

 

  运行结果:

  嗯,和上面MATLAB还是一样的,看来代码没问题咯。但是呢,从输出可以看到,这里一次拟合只能插一个值,那么如果同样的数据,我们想要插入两个点,这个代码就行不通了。其实改也很简单,只要挪一下do-while循环的位置就可以了。下面是改进的代码。

#include<iostream>
#include<string>
#include<vector>
using namespace std;

// 函数声明,是的其在完整定义前可以被(main函数)使用
double Lagrange(int N,vector<double>&X,vector<double>&Y,double x);

int main()
{
	char a='n';
	
	cout<<"请输入插值次数n的值:"<<endl;
	int N;
	cin>>N;
	vector<double>X(N,0);
	vector<double>Y(N,0);
	cout<<"请输入插值点对应的值及其函数值(Xi,Yi):"<<endl;
	for(int a=0;a<N;a++){
		cin>>X[a]>>Y[a];}
	// 注意do从上面移动到了这里,也就是改变了循环的内容
	do{
		cout<<"请输入要求值的x的值:"<<endl;
		double x;
		cin>>x;
		double result=Lagrange(N,X,Y,x);
		cout<<"由拉格朗日插值法得出的结果:"<<result<<endl;
		cout<<"是否要继续?<y/n>:";
		cin>>a;
	}while(a=='y');
	return 0;
}

//vector<double>,创建空容器,其对象类型为double类
double Lagrange(int N,vector<double>&X,vector<double>&Y,double x)
{
	double result=0;
	for(int i=0;i<N;i++){
		double temp=Y[i];
		for(int j=0;j<N;j++){
			if(i != j){
				temp = temp*(x-X[j]);
				temp = temp/(X[i]-X[j]);
			}
		}// 与上文MATLAB实现不同,这里直接将拉格朗日多项式的各项直接相加
		result += temp;
	}
	return result; 
}

  运行结果: