三次样条插值——三弯矩法（C语言实现）

1、原理

2、案例

3、代码实现

#include<iostream>
#include<iomanip>
using namespace std;
#define max 50 

float x[max], y[max], h[max];//Define a specific array to store the original data
float c[max], a[max], fm[max];
float f(int x1, int x2, int x3)     
{
	float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);   //Catch up method
	float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);
	return (a - b) / (x[x3] - x[x1]);
}   
int main()     //Find difference
{
	void print(int n);    //Function declaration
	void chaf(int n);     //Function declaration
	int n, i;               
	char ch;
	cout << "Enter the maximum subscript of x:";
	cin >> n;
	for (i = 0; i <= n; i++)
	{
		cout << "Please enterx[" << i << "]=";
		cin >> x[i];
	}
	for (i = 0; i <= n; i++)
	{
		cout << "Please entery[" << i << "]=";
		cin >> y[i];	
	}
	cout << endl;
	do
	{
		for (i = 0; i < n; i++)            //Find step size
			h[i] = x[i + 1] - x[i];
		cout << "Please enter boundary conditions:";
		int p;
		float f0, f1;
		cin >> p;
		switch (p)
		{
		case 1:
			cout << "inputY0\' Y" << n << "\'\n";
			cin >> f0 >> f1;
			c[0] = 1;
			a[n] = 1;
			fm[0] = 6 * ((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];
			fm[n] = 6 * (f1 - (y[n] - y[n - 1]) / (x[n] - x[n - 1])) / h[n - 1];
			break;
		case 2:
			cout << "inputY0\" Y" << n << "\"\n";
			cin >> f0 >> f1;
			c[0] = a[n] = 0;
			fm[0] = 2 * f0; fm[n] = 2 * f1;
			break;
		default:
			cout << "Not available\n";//undetermined
		}
		for (i = 1; i < n; i++)
			fm[i] = 6 * f(i - 1, i, i + 1);
		for (i = 1; i < n; i++)
		{
			a[i] = h[i - 1] / (h[i] + h[i - 1]);
			c[i] = 1 - a[i];
		}
		a[n] = h[n - 1] / (h[n - 1] + h[n]);
		chaf(n);            //Call function
		cout << "Interpolation function：";
		print(n);           //Call function
		cin >> ch;
	} 
while (ch == 'y' || ch == 'Y');
system("pause");
return 0;
}
void print(int n)
{
	cout << setprecision(6);
	for (int i = 0; i < n; i++)
	{
		cout << i + 1 << ": (" << x[i] << " , " << x[i + 1] << ")\n" << "\t";
		cout << "S" << i + 1 << "=";
		float t = fm[i] / (6 * h[i]);
		if (t > 0)
			cout << -t << "*(x-" << x[i + 1] << ")^3";
		else
			cout << -t << "*(x-" << x[i + 1] << ")^3";
		t = fm[i + 1] / (6 * h[i]);
		if (t > 0)
			cout << " + " << t << "*(x-" << x[i] << ")^3";
		else
			cout << "-" << t << "*(x-" << x[i] << ")^3";
		cout << "\n\t";
		t = (y[i] - fm[i] * h[i] * h[i] / 6) / h[i];
		if (t > 0)
			cout << "- " << t << "*(x-" << x[i + 1] << ")";
		else
			cout << "- " << -t << "*(x-" << x[i + 1] << ")";
		t = (y[i + 1] - fm[i + 1] * h[i] * h[i] / 6) / h[i];
		if (t > 0)
			cout << " + " << t << "*(x-" << x[i] << ")";
		else
			cout << "-" << -t << "*(x-" << x[i] << ")";
		cout << endl;
	}
	cout << endl;
}
void chaf(int n)
{
	float B[max];
	B[0] = c[0] / 2;
	int i;
	for (i = 1; i < n; i++)
		B[i] = c[i] / (2 - a[i] * B[i - 1]);
	fm[0] = fm[0] / 2;
	for (i = 1; i <= n; i++)
		fm[i] = (fm[i] - a[i] * fm[i - 1]) / (2 - a[i] * B[i - 1]);
	for (i = n - 1; i >= 0; i--)
		fm[i] = fm[i] - B[i] * fm[i + 1];
}

4、结果

 5、误差分析与心得体会