《曲线拟合 给定10个点 拟合曲线 复制即用 c++》由会员分享，可在线阅读，更多相关《曲线拟合 给定10个点 拟合曲线 复制即用 c++（11页珍藏版）》请在金锄头文库上搜索。
1、#include #include #include #include#defineF1(t)(1+(t)*(t)*(2*(t)-3)#defineF2(t)(t)*(t)*(-2*(t)+3)#defineF3(t)(t)*(1+(t)*(t)-2)#defineF4(t)(t)*(t)*(t)-1)using namespace std;void tangent(double *y,double *dy,int n) /计算各个节点的切矢，采用自由端点条件 /追赶法求解三对角方程组int i;double *a,*b,*c,*d;a=new doublen;b=new doublen;c=
2、new doublen;d=new doublen;/追赶法计算系数for(i=1;i n;i+)ai=1;b0=2;bn-1=2;for(i=1;i n-1;i+)bi=4;for(i=0;i n-1;i+)ci=1;d0=3*(y1-y0); dn-1=3*(yn-1-yn-2);for(i=1;i n-1;i+)di=3*(yi+1-yi-1);/追赶法回代 double *L,*U; L=new doublen; U=new doublen; U0=b0;for(i=1;i n;i+)Li=ai/Ui-1;Ui=bi-Li*ci-1;double *yy; yy=new doublen
3、;yy0=d0;for(i=1;i =0;i-)dyi=(yyi-ci*dyi+1)/Ui;/delete a,b,c,d,L,U,yy;xa,int n)double spline(double *x,double *y,double *dy,double/三次 Hermite 插值数据点 xn,yn,切矢 dyn返回xa处的插值数值要求x严格递增，刈0 =xa =xn-1int i;for(i=0;i n-1;i+)if(xa =xi+1)break;double t = (xa-xi)/(xi+1-xi);return F1(t)*yi+F2(t)*yi+1+F3(t)*dyi+F4(t
4、)*dyi+1;int main()doubles_ua,xa_ua;/A 类型循环变量doubles_ub,xa_ub;/B 类型循环变量doubles_uc,xa_uc; C类型循环变量doubles_a_result,xa_a_result;/doubles_b_result,xa_b_result;/doubleaverage_s,average_xa;doubledy_ua4;doubledy_ub4;doubledy_uc4;doublex_ua4=0.000,0,0,0;doubley_ua4=0.0000,0,0,0;doublex_ub4=0.000,0,0,0;double
5、y_ub4=0.0000,0,0,0;double x_uc4=0.000,0,0,0;double y_uc4=0.0000,0,0,0;/DNA/double X100.000,7.000,11.197,17.910,28.649,45.826,73.301,117.251,187.551,300.00;/double Y10 = 272,3008,5854,10287,18398,32204,56751,95727,151346,232604; /HCG/double X100.000,0.500,1.000,5.000,10.000,50.000,100.000,500.000,100
6、0.000,5000.000;/double Y10 = 1477,1869,1920,4279,7010,28345,55156,245217,415291,1294682; /HSV-M/double X10 = 0.000,0.500,0.645,0.832,1.073,1.385,1.786,2.304,2.973,3.835;/double Y10 = 349,4273,5487,7054,9072,11685,15043,19381,24984,32203;/Trepdouble X10 = 0.000,0.400,0.761,1.450,2.760,5.254,10.001,19
7、.039,36.245,70.000; double Y10 = 988,2620,4889,9686,19423,39364,81962,152041,249416,398202; int n;/*string s_a;用于 a b c 段线string s_b;用于 b c 段线string s_c;用于c段线string s_avarage;用于b c段线 求均值*/for(n=0;n9;n+)/computing(X,Y,n);int k_a,k_b,k_c;/插值函数调用，并把结果输入到文件if(n = 0)coutnendl;for(k_a=0;k_a4;k_a+)x_uak_a
8、= Xk_a+n;y_uak_a = Yk_a+n;/给x、y赋初值。ofstream ofs_all(d:spline_all.txt);写文件流if(!ofs_all) return 10000; tangent(y_ua,dy_ua,4);for(xa_ua = X0;xa_ua = X1;xa_ua += 0.01)s_ua=spline(x_ua,y_ua,dy_ua,xa_ua,4);/fprintf(ofs_all,%f#%fn ,xa,s); ofs_allsetw(15)xa_uasetw(15)s_uaendl;/ofs_all.close();/写入第一部分else if
9、(n = 1)coutnendl;for(k_a=0;k_a4;k_a+)x_uak_a = Xk_a;y_uak_a = Yk_a;/给x、y赋初值。tangent(y_ua,dy_ua,4);ofstream ofs_all(d:spline_all.txt,ios_base:app);写文件流if(!ofs_all) return 10000;for(xa_ua = X1;xa_ua = X2;xa_ua += 0.01)s_ua=spline(x_ua,y_ua,dy_ua,xa_ua,4);xa_a_result = xa_ua;s_a_result= s_ua;for(k_b=0;
10、k_b4;k_b+)x_ubk_b = Xk_b+1;y_ubk_b = Yk_b+1;tangent(y_ub,dy_ub,4);/for(xa_ub = X1;xa_ub = X2;xa_ub += 0.01)/xa_ub = xa_ua-0.01;s_ub=spline(x_ub,y_ub,dy_ub,xa_ua,4);average_s = (s_a_result+s_ub)/2;average_xa = (xa_a_result+xa_ub)/2;ofs_allsetw(15)average_xasetw(15)average_sendl;/ofs_all.close();/C x2
11、x3else if(n = 2)coutnendl;for(k_a=0;k_a4;k_a+)x_uak_a = Xk_a+n-2;y_uak_a = Yk_a+n-2;/给x、y赋初值。tangent(y_ua,dy_ua,4);ofstream ofs_all(d:spline_all.txt,ios_base:app);写文件流if(!ofs_all) return 10000;for(xa_ua = X2;xa_ua = X3;xa_ua += 0.01)s_ua=spline(x_ua,y_ua,dy_ua,xa_ua,4); xa_a_result = xa_ua;s_a_resul
12、t = s_ua;/bfor(k_b=0;k_b4;k_b+)x_ubk_b = Xk_b+n-1;y_ubk_b = Yk_b+n-1; tangent(y_ub,dy_ub,4);xa_ub = xa_ua-0.01;xa_uc = xa_ub;s_ub=spline(x_ub,y_ub,dy_ub,xa_ua,4);/average_s = (s_a_result+s_ub)/2;/average_xa = (xa_a_result+xa_ub)/2;/cfor(k_c=0;k_c4;k_c+)x_uck_c = Xk_c+n;y_uck_c = Yk_c+n;tangent(y_uc,
13、dy_uc,4); s_uc=spline(x_uc,y_uc,dy_uc,xa_ua,4); average_s = (s_a_result+s_ub+s_uc)/3; average_xa = (xa_a_result+xa_ub+xa_uc)/3;ofs_allsetw(15)average_xasetw(15)average_sendl;/ ofs_all.close();else if(n = 3)coutnendl;for(k_a=0;k_a4;k_a+)x_uak_a = Xk_a+n-2;y_uak_a = Yk_a+n-2;/给x、y赋初值。tangent(y_ua,dy_ua,4);ofstream ofs_all(d:spline_all.txt,ios_base:app);写文件流 if(!ofs_all) return 10000;for(xa_ua = X3;xa_ua = X4;xa_u
