/* 测试数据来自课本55页例5 《数值分析》清华大学出版社第四版*/
//输入
3
27.7 4.1
28 4.3
29 4.1
30 3.0
1
3.0 -4.0
//输出
输出三次样条插值函数：
1: [27.7 , 28]
13.07*(x - 28)^3 + 0.22*(x - 27.7)^3
+ 14.84*(28 - x) + 14.31*(x - 27.7)
2: [28 , 29]
0.066*(29 - x)^3 + 0.1383*(x - 28)^3
+ 4.234*(29 - x) + 3.962*(x - 28)
3: [29 , 30]
0.1383*(30 - x)^3 - 1.519*(x - 29)^3
+ 3.962*(30 - x) + 4.519*(x - 29)
//三次样条插值函数
#include<iostream>
#include<iomanip>
using namespace std;
const int MAX = 50;
float x[MAX], y[MAX], h[MAX];
float c[MAX], a[MAX], fxym[MAX];
float f(int x1, int x2, int x3){
float a = (y[x3] - y[x2]) / (x[x3] - x[x2]);
float b = (y[x2] - y[x1]) / (x[x2] - x[x1]);
return (a - b)/(x[x3] - x[x1]);
} //求差分
void cal_m(int n){ //用追赶法求解出弯矩向量M……
float B[MAX];
B[0] = c[0] / 2;
for(int i = 1; i < n; i++)
B[i] = c[i] / (2 - a[i]*B[i-1]);
fxym[0] = fxym[0] / 2;
for(i = 1; i <= n; i++)
fxym[i] = (fxym[i] - a[i]*fxym[i-1]) / (2 - a[i]*B[i-1]);
for(i = n-1; i >= 0; i--)
fxym[i] = fxym[i] - B[i]*fxym[i+1];
}
void printout(int n);
int main(){
int n,i; char ch;
do{
cout<<"Please put in the number of the dots:";
cin>>n;
for(i = 0; i <= n; i++){
cout<<"Please put in X"<<i<<':';
cin>>x[i]; //cout<<endl;
cout<<"Please put in Y"<<i<<':';
cin>>y[i]; //cout<<endl;
}
for(i = 0; i < n; i++) //求 步长
h[i] = x[i+1] - x[i];
cout<<"Please 输入边界条件\n 1: 已知两端的一阶导数\n 2:两端的二阶导数已知\n 默认:自然边界条件\n";
int t;
float f0, f1;
cin>>t;
switch(t){
case 1:cout<<"Please put in Y0\' Y"<<n<<"\'\n";
cin>>f0>>f1;
c[0] = 1; a[n] = 1;
fxym[0] = 6*((y[1] - y[0]) / (x[1] - x[0]) - f0) / h[0];
fxym[n] = 6*(f1 - (y[n] - y[n-1]) / (x[n] - x[n-1])) / h[n-1];
break;
case 2:cout<<"Please put in Y0\" Y"<<n<<"\"\n";
cin>>f0>>f1;
c[0] = a[n] = 0;
fxym[0] = 2*f0; fxym[n] = 2*f1;
break;
default:cout<<"不可用\n";//待定
};//switch
for(i = 1; i < n; i++)
fxym[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]);
cal_m(n);
cout<<"\n输出三次样条插值函数：\n";
printout(n);
cout<<"Do you to have anther try ? y/n :";
cin>>ch;
}while(ch == 'y' || ch == 'Y');
return 0;
}
void printout(int n){
cout<<setprecision(6);
for(int i = 0; i < n; i++){
cout<<i+1<<": ["<<x[i]<<" , "<<x[i+1]<<"]\n"<<"\t";
/*
cout<<fxym[i]/(6*h[i])<<" * ("<<x[i+1]<<" - x)^3 + "<<<<" * (x - "<<x[i]<<")^3 + "
<<(y[i] - fxym[i]*h[i]*h[i]/6)/h[i]<<" * ("<<x[i+1]<<" - x) + "
<<(y[i+1] - fxym[i+1]*h[i]*h[i]/6)/h[i]<<"(x - "<<x[i]<<")\n";
cout<<endl;*/
float t = fxym[i]/(6*h[i]);
if(t > 0)cout<<t<<"*("<<x[i+1]<<" - x)^3";
else cout<<-t<<"*(x - "<<x[i+1]<<")^3";
t = fxym[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] - fxym[i]*h[i]*h[i]/6)/h[i];
if(t > 0)cout<<"+ "<<t<<"*("<<x[i+1]<<" - x)";
else cout<<"- "<<-t<<"*("<<x[i+1]<<" - x)";
t = (y[i+1] - fxym[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<<endl;
}
cout<<endl;
}
评论0