zuoye.cpp

用C++实行数值插值运算并用matlab绘图

#include <iostream.h>
#include <fstream.h>
#include <cmath >
using namespace std;
int main()
{const int nodes=11;//定义符号常量nodes，大小与节点数相等
const int insertnodes=30;//定义符号常量insernodes，大小与插值点数相等

	double xx[nodes]={-1,-0.8,-0.6,-0.4,-0.2,0,0.2,0.4,0.6,0.8,1};//定义放节点的数组xx[ ]
	double yy[nodes]={1.96E-02,3.03E-02,5.26E-02,1.111E-01,3.333E-01,1.00E-00,3.333E-01,1.111E-01,5.26E-02,3.03E-02,1.96E-02
};//定义放节点的数组yy[ ]
double   
xi[insertnodes]={-0.95,-0.9,-0.85,-0.75,-0.7,-0.65,-0.55,-0.5,-0.45,-0.35,-0.3,-0.25, -0.15,-0.1,-0.05,0.05,0.1,0.15,0.25,0.3,0.35,0.45,0.5,0.55,0.65,0.7,0.75,0.85,0.9,0.95};//设置插值点的x值
double x[nodes+insertnodes];    double l,L;    int i, k,m,n;  int j=0;    double f;
	//通过两个for循环把所有节点和插值点放到数组x[ ]
	
	for (i=0;i<nodes;i++)
	{
		x[j]=xx[i];
		j++;
	}
	for (i=0;i<insertnodes;i++)
	{
		x[j]=xi[i];
		j++;}
		//把数组x[ ]里的数从大到小排列
	for (i=1;i<nodes+insertnodes-1;i++)
	{
		for (j=0;j<nodes+insertnodes-i;j++)
		{
			if (x[j]<x[j+1])
			{
				l=x[j+1];
				x[j+1]=x[j];
				x[j]=l;
			}
		}
	}

	
	ofstream outfile("E://zuo.m");//输出个zuo.m文件到E://zuo.m
	//输出zuo.m文件中的x数组
	outfile<<"x=[";
	for (k=0;k<nodes+insertnodes-1;k++)
		outfile<<x[k]<<",";
	outfile<<x[k]<<"];"<<endl;
	//输出zuo.m文件中的y数组
	outfile<<"y=[";
	//对数组中的所有x值进行处理,如果x值是节点,直接取它的y值,是插值点,用拉格朗日插值法求出其y值
	n=0;//用于控制节点数组的下标
	for (m=0;m<nodes+insertnodes;m++)//对数组x循环操作(数组x包含了节点和插值点的x值)
	{
		if (x[m]==xx[n])//说明是节点,直接取其y值
		{
			L=yy[n];//取节点的y值放到L中
			//cout<<"节点:"<<x[m]<<"值:"<<L<<endl;
			n++;//指向下一个节点的x值
		}
		else//说明是插值点
		{
			//用拉格朗日插值法求出插值点的y值
			L=0;
			for (i=0;i<nodes;i++)
			{
				l=1.0;
				for (j=0;j<nodes;j++)
					if (i!=j) l=l*(x[m]-xx[j])/(xx[i]-xx[j]);
				L=L+l*yy[i];
			}
			//显示在屏幕上便于对比查看计算结果
			cout<<"插值点:"<<x[m]<<"\t";
			cout<<"拉氏值:"<<L<<"\t";//
			f=1.0/(1+50*pow(x[m],2));//理论函数的精确值
			cout<<"理论值:"<<f<<"\n";
		}
		//把y值输出到zuo.m文件中
		if (m<nodes+insertnodes-1)
			outfile<<L<<",";
		else
			outfile<<L<<"];"<<endl;
	}
	outfile<<"x1=-1:0.01:1;"<<endl;//用MatLab绘制出y=1/(1+50*x^2)的理论曲线
	outfile<<"y1=1./(1+50*x1.^2);"<<endl;
	outfile<<"plot(x1,y1,x,y,':');"<<endl;//wei.m文件绘图(同时绘出理论曲线和插值曲线)
	outfile<<"legend('理论曲线y=1/(1+50*x^2)','拉氏插值曲线');"<<endl;
	outfile.close();//关闭zuo.m文件
	return 0;
}