实例介绍
【实例截图】
【核心代码】
// 三次样条曲线拟合.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; }
标签: 计算
小贴士
感谢您为本站写下的评论,您的评论对其它用户来说具有重要的参考价值,所以请认真填写。
- 类似“顶”、“沙发”之类没有营养的文字,对勤劳贡献的楼主来说是令人沮丧的反馈信息。
- 相信您也不想看到一排文字/表情墙,所以请不要反馈意义不大的重复字符,也请尽量不要纯表情的回复。
- 提问之前请再仔细看一遍楼主的说明,或许是您遗漏了。
- 请勿到处挖坑绊人、招贴广告。既占空间让人厌烦,又没人会搭理,于人于己都无利。
关于好例子网
本站旨在为广大IT学习爱好者提供一个非营利性互相学习交流分享平台。本站所有资源都可以被免费获取学习研究。本站资源来自网友分享,对搜索内容的合法性不具有预见性、识别性、控制性,仅供学习研究,请务必在下载后24小时内给予删除,不得用于其他任何用途,否则后果自负。基于互联网的特殊性,平台无法对用户传输的作品、信息、内容的权属或合法性、安全性、合规性、真实性、科学性、完整权、有效性等进行实质审查;无论平台是否已进行审查,用户均应自行承担因其传输的作品、信息、内容而可能或已经产生的侵权或权属纠纷等法律责任。本站所有资源不代表本站的观点或立场,基于网友分享,根据中国法律《信息网络传播权保护条例》第二十二与二十三条之规定,若资源存在侵权或相关问题请联系本站客服人员,点此联系我们。关于更多版权及免责申明参见 版权及免责申明
网友评论
我要评论