sanciyangtiaochazhi.rar_sanciyangtiaochazhi_一阶导数_三次样条插值_三次自然样条_样

/* 测试数据来自课本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; }