三次样条曲线拟合 三次样条曲线拟合

// 三次样条曲线拟合.cpp : Defines the entry point for the console application. // #include "stdafx.h" #include<iostream.h> #include<iomanip.h> #define N 50 int main(int argc, char* argv[]) { int i,j,k,n; double X[N],Y[N],h[N],u[N],d[N],m[N],a[N],b[N],A[N][N],S[N][4]; void line_equations(double A[][50],double X[],int n=1); //变量的定义 说明: k=0 to n-1 u[0]+=m[k]*a[k]; // k=0 to n-1 u[n-1]+=m[k]*b[k]; 为附加的两个条件 cout<<"请输入点的个数n:"; cin>>n; cout<<"请输入各个点的坐标:"<<endl; for(i=0;i<=n-1;i++){ cout<<"X["<<i<<"]="; cin>>X[i]; cout<<"Y["<<i<<"]="; cin>>Y[i]; }; //坐标的输入 cout<<"请输入两个附加条件:m[k]的线性组合"<<endl; for(i=0;i<=n-1;i++){ cout<<"+m["<<i<<"]*"; cin>>a[i]; }; cout<<"="; cin>>u[0]; for(i=0;i<=n-1;i++){ cout<<"+m["<<i<<"]*"; cin>>b[i]; }; cout<<"="; cin>>u[n-1]; //附加条件的输入 for(k=0;k<=n-2;k++){ h[k]=X[k+1]-X[k]; d[k]=(Y[k+1]-Y[k])/(X[k+1]-X[k]); }; if(n>=3){ for(k=1;k<=n-2;k++){ u[k]=6*(d[k]-d[k-1]); }; for(i=1;i<=n-2;i++){ for(j=0;j<=n;j++){ A[i][j]=0; } }; }; //初始化方程组的系数增广矩阵 for(i=0;i<=n-1;i++){ A[0][i]=a[i]; A[n-1][i]=b[i]; A[i][n]=u[i]; }; if(n>=3){ for(i=1;i<=n-2;i++){ A[i][i-1]=h[i-1]; A[i][i]=2*(h[i-1]+h[i]); A[i][i+1]=h[i]; } }; //方程组的系数增广矩阵的填充 line_equations(A,m,n); //求解m[k]; for(k=0;k<=n-2;k++){ S[k][0]=m[k]*X[k+1]*X[k+1]*X[k+1]/(6*h[k])+(Y[k]/h[k]-m[k]*h[k]/6)*X[k+1]; S[k][0]-=m[k+1]/(6*h[k])*X[k]*X[k]*X[k]+(Y[k+1]/h[k]-m[k+1]*h[k]/6)*X[k]; S[k][1]=m[k+1]*X[k]*X[k]/(2*h[k])+(Y[k+1]/h[k]-m[k+1]*h[k]/6); S[k][1]-=m[k]*X[k+1]*X[k+1]/(2*h[k])+(Y[k]/h[k]-m[k]*h[k]/6); S[k][2]=m[k]*X[k+1]/(2*h[k])-m[k+1]*X[k]/(2*h[k]); S[k][3]=(m[k+1]-m[k])/(6*h[k]); }; //求解三次样条曲线的各次方的系数 cout<<"下面是三次样条曲线:"<<endl; for(i=0;i<=n-2;i++){ cout<<"S["<<i<<"]["<<"X]= "<<S[i][3]<<" *X^3 "; if(S[i][2]>=0) cout<<"+"; cout<<S[i][2]<<" *X^2 "; if(S[i][1]>=0) cout<<"+"; cout<<S[i][1]<<" *X "; if(S[i][0]>=0) cout<<"+"; cout<<S[i][0]<<endl; }; //输出三次样条曲线多项式 cin.get(); cin.get(); return 0; }; void line_equations(double A[][50],double X[],int n=1){ int i,j,k,s=0; double temp; if(n==1 && A[0][0]!=0){ X[0]=A[0][1]/A[0][0]; return ; }; //为一阶时直接求出解 for(k=0;k<=n-2;k++){ for(i=k;i<=n;i++){ if(A[i][k]!=0) break; }; //求出从上往下第一个不为0的行 if(i>=n){ cout<<"无解!"; return ; }; //若在化归上三角矩阵的时候碰到有缺阶数的时候，此时无解或有无穷多解，统一为无解 if(i!=k){ for(j=k;j<=n;j++){ temp=A[k][j]; A[k][j]=A[i][j]; A[i][j]=temp; }; }; //第一个不为0的行来与k行交换，以免除数为0 for(i=k+1;i<=n-1;i++){ for(j=k+1;j<=n;j++){ A[i][j]-=((A[i][k]/A[k][k])*A[k][j]); }; A[i][k]=0; }; //第k次循环时，将第1列化为0 }; //至此，将方程组的系数增广矩阵化为了上三角阵 X[n-1]=A[n-1][n]/A[n-1][n-1]; //从下向上迭代求出X[k],先求最后一个 for(i=n-2;i>=0;i--){ for(j=n-2;j>=i;j--){ A[i][n]-=X[j+1]*A[i][j+1]; }; X[i]=A[i][n]/A[i][i]; }; //迭代求出X[k] return; }