three.cpp

计算机数值分析实验题目 三次样条插值的程序 满足三次样条插值的三弯矩法 编制第一与第二边界条件的程序

#include <iostream>//数据流输入／输出
#include <fstream> //文件输入／输出
#include <iomanip>//参数化输入／输出
#include <cmath>  /定义数学函数


using namespace std;

struct Point
{
	double x, y;
};

const int MaxPoint = 50;

void PrintMatrix(double Matrix[MaxPoint][MaxPoint + 1], int NumOfLine)
{
	int ci, cj;
	for(ci = 0; ci < NumOfLine; ci++)
	{
		for(cj = 0; cj < NumOfLine + 1; cj++)
			cout << setprecision(5) <<  setw(10) << Matrix[ci][cj];
		cout << endl;
	}
	return;
}

// 解矩阵 返回0表示正常 1表示出错 高斯消元法
int WorkOutMatrix(double Matrix[MaxPoint][MaxPoint + 1], int NumOfLine)
{
	int ci, cj, cm;  // 循环变量
	double temp;
	// 化简成上三角矩阵
	for(ci = 0; ci < NumOfLine; ci++)
	{
		// 找到一个不是0的行进行交换
		for(cj = ci; cj < NumOfLine; cj++)
			if(abs(Matrix[cj][ci]) > 1e-10)
				break;
		if(cj == NumOfLine)
			continue;
		// 交换 ci cj 行
		if(ci != cj)
		{
			for(cm = 0; cm <= NumOfLine; cm++)
			{
				temp = Matrix[ci][cm];
				Matrix[ci][cm] = Matrix[cj][cm];
				Matrix[cj][cm] = temp;
			}
		}
		// 消元
		for(cj = ci + 1; cj < NumOfLine; cj++)
		{
			if(abs(Matrix[cj][ci]) < 1e-10)
				continue;
			temp = 0.0 - Matrix[cj][ci] / Matrix[ci][ci];
			// 把第ci行乘加到cj行
			for(cm = ci; cm <= NumOfLine; cm++)
			{
				//cout << "## " << Matrix[cj][cm] << " += " << Matrix[ci][cm] << "* -" << Matrix[cj][ci] << " / " << Matrix[ci][ci] 
				//<< "  temp = " << temp << endl; 
				Matrix[cj][cm] += Matrix[ci][cm] * temp;
				if(abs(Matrix[cj][cm]) < 1e-10)
					Matrix[cj][cm] = 0;
			}
		}
		//cout << "=======================" << endl;
		//PrintMatrix(Matrix, NumOfLine);
		//cout << "=======================" << endl;

	}
	// 化成单位阵
	for(ci = NumOfLine - 1; ci >= 0; ci--)
	{
		if(Matrix[ci][ci] < 1e-4)
			return 1;
		for(cj = ci - 1; cj >= 0; cj--)
		{
			temp = 0 - Matrix[cj][ci] / Matrix[ci][ci];
			Matrix[cj][ci] += Matrix[ci][ci] *temp;
			Matrix[cj][NumOfLine] += Matrix[ci][NumOfLine] * temp;
			if(abs(Matrix[cj][ci]) < 1e-4)
				Matrix[cj][ci] = 0;
			if(abs(Matrix[cj][NumOfLine]) < 1e-4)
				Matrix[cj][NumOfLine] = 0;
		}
	}
	for(ci = 0; ci < NumOfLine; ci++)
	{
		Matrix[ci][NumOfLine] /= Matrix[ci][ci];
		Matrix[ci][ci] = 1;
	}
	return 0;
}
int main()
{
	int ci;
	int NumOfPoint;
	double der1,der2;  // 边界条件
	Point p[MaxPoint];
	double Matrix[MaxPoint][MaxPoint + 1] = {0};
	double u[MaxPoint] = {0};
	double l[MaxPoint] = {0};
	double d[MaxPoint] = {0};
	ifstream in("in.txt");

	// 数据输入
	//cout << "三次样条插值" << endl;
	//cout << "============" << endl;
	//cout << "输入坐标点数目:";
	//cin >> NumOfPoint;
	in >> NumOfPoint;

	for(ci = 0; ci < NumOfPoint; ci++)
	{
		//cout << "输入第" << ci + 1 << "个点x坐标:" << endl;
		//cin >> p[ci].x;
		//cout << "输入第" << ci + 1 << "个点y坐标:" << endl;
		//cin >> p[ci].y;
		in >> p[ci].x;
		in >> p[ci].y;
	}
	//cout << "输入边界条件S'(" << p[0].x << ") = ";
	//cin >> der1;
	//cout << "输入边界条件S'(" << p[NumOfPoint - 1].x << ") = ";
	//cin >> der2;
	in >> der1;
	in >> der2;

	// 数据写入矩阵
	for(ci = 0; ci < NumOfPoint; ci++)  // 写入对角线
		Matrix[ci][ci] = 2;
	
	for(ci = 1; ci < NumOfPoint - 1; ci++)  // 写入另外两个对角线
	{
		u[ci] = (p[ci].x - p[ci - 1].x) / (p[ci + 1].x - p[ci - 1].x);
		l[ci] = (p[ci + 1].x - p[ci].x) / (p[ci + 1].x - p[ci - 1].x);
	}
	u[NumOfPoint - 1] = 1;
	l[0] = 1;
	for(ci = 1; ci < NumOfPoint - 1; ci++)
	{
		Matrix[ci][ci - 1] = u[ci];
		Matrix[ci][ci + 1] = l[ci];
	}
	Matrix[0][1] = l[0];
	Matrix[NumOfPoint - 1][NumOfPoint - 2] = u[NumOfPoint - 1];

	// 将d的值写入矩阵
	for(ci = 1; ci < NumOfPoint - 1; ci++)
		d[ci] = 6 * ((p[ci + 1].y - p[ci].y) / (p[ci + 1].x 
			- p[ci].x) - (p[ci].y - p[ci - 1].y) / (p[ci].x - p[ci - 1].x))
			/ (p[ci + 1].x - p[ci - 1].x);
	d[0] = 6 * ((p[1].y - p[0].y) / (p[1].x - p[0].x) - der1) / (p[1].x - p[0].x);
	d[NumOfPoint - 1] = 6 * (der2 - (p[NumOfPoint - 1].y - p[NumOfPoint - 2].y) 
		/ (p[NumOfPoint - 1].x - p[NumOfPoint - 2].x)) / (p[NumOfPoint - 1].x - p[NumOfPoint - 2].x);
	for(ci = 0; ci < NumOfPoint; ci++)
		Matrix[ci][NumOfPoint] = d[ci];

	// 测试
	//PrintMatrix(Matrix, NumOfPoint);
	WorkOutMatrix(Matrix, NumOfPoint);
	//cout << "=======================" << endl;
	//PrintMatrix(Matrix, NumOfPoint);
	cout.setf(ios::left);
	for(ci = 0; ci < NumOfPoint; ci++)
		cout << "M" << ci << " = " << setprecision(5) <<  setw(10) << Matrix[ci][NumOfPoint] << '\t';
	cout << endl;

	return 0;
}