第二型样条插值.cpp

第一型样条插值

/*
	第一型样条插值，边界条件为两端的二阶导数值已知


*/
#include <iostream>
#include <vector>
#include <fstream>
#include "zhuigang.h"
#include "conio.h"//_getch()

using namespace std;

void main()
{
	vector<double> x, y;
	int n;//数值的个数
	fstream filex("Spline2x.txt");//数据文件
	fstream filey("Spline2y.txt");
	double temp;
	while (filex >> temp)
		x.push_back(temp);
	while (filey >> temp)
		y.push_back(temp);

	n = x.end() - x.begin();
	double df0, dfn;
	df0 = y[n];//端点的导数值
	dfn = y[n + 1];

	vector<double> h, mui, b, lmel, M, d;

	for(int i = 0; i < n - 1; i++)//数据间隔
	{
		h.push_back(x[i + 1] - x[i]);
	}

	for(i = 0; i < n - 2; i++)//μ值
	{
		mui.push_back(h[i] / (h[i] + h[i + 1]));
	}
	mui.push_back(1);

	for(i = 0; i < n - 2; i++)//λ值
	{
		lmel.push_back(1 - mui[i]);
	}

	for(i = 0; i < n - 2; i++)//d值
	{
		temp = y[i] / h[i] / (h[i] + h[i + 1]) - y[i + 1] / h[i] / h[i + 1]
			+ y[i + 2] / (h[i] + h[i + 1]) / h[i + 1];
		d.push_back(6 * temp);
	}
	d[0] -= mui[0] * df0;
	d[n - 3] -= lmel[n - 3] * dfn;

	for(i = 0; i < n - 2; i++)
	{
		b.push_back(2);
	}
	
	vector<double> a, c, Mlast;//Mlast第一和最后一个数值已知，其它部分由M组成
	for(i = 0; i < n - 3; i++)
	{
		a.push_back(mui[i + 1]);
		c.push_back(lmel[i]);
	}
	M = zhuigang(a, b, c, d);//追赶法求M值

	Mlast.push_back(df0);
	for(i = 0; i < n - 2; i++)
	{
		Mlast.push_back(M[i]);
	}
	Mlast.push_back(dfn);

	for (i = 0; i < n; i++)
	{
		cout << Mlast[i] << "\t";
	}

	do//检查程序输出
	{
		cout << "\ninput x:";
		cin >> temp;

		int j = 0;//找出x所在的区间
		for(i = 0; i < n, j == 0; i++)
		{
			if(x[i] > temp)
				j = i + 1;
		}
		j -= 2;
		if(j < 0 || j > n - 1) 
		{
			cout << "the input is out of range";
			goto another;
		}

		cout << "S(" << temp << ") is: ";

		temp =y[j] + ((y[j + 1] - y[j]) / h[j] - (Mlast[j] / 3 + Mlast[j + 1] / 6) * h[j]) 
			 * (temp - x[j]) + Mlast[j] / 2 * (temp - x[j]) * (temp - x[j])
			+ (Mlast[j + 1] - Mlast[j]) / 6 / h[j]
			 * (temp - x[j]) * (temp - x[j]) * (temp - x[j]);

		cout << temp << endl;
another:cout << "press any key to see another S(x) except x(exit)";

	}while(_getch() != 'x');

}