数值计算课程设计

1、1、 经典四阶龙格库塔法解一阶微分方程1.1、 算法说明龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。4阶龙格-库塔方法(RK4)可模拟N=4的泰勒方法的精度。这种算法可以描述为，自初始点开始，利用下面的计算方法生成近似序列(1-1)1.2、 经典四阶龙格库塔法解一阶微分方程算法流程图图1-1 经典四阶龙格库塔法解一阶微分方程算法流程图- 1 - / 611.3、 经典四阶龙格库塔法解一阶微分方程程序调试图1-2 经典四阶龙格库塔法解一阶微分方程程序调试1.4

2、、 经典四阶龙格库塔法解一阶微分方程代码#include <iostream>#include <iomanip>using namespace std;/f为函数的入口地址，x0、y0为初值，xn为所求点，step为计算次数double Runge_Kuta( double (*f)(double x, double y), double x0, double y0, double xn, int step ) double k1,k2,k3,k4,result; double h=(xn-x0)/step; if(step<=0) return(y0); if

3、(step=1) k1=f(x0,y0); k2=f(x0+h/2, y0+h*k1/2); k3=f(x0+h/2, y0+h*k2/2); k4=f(x0+h, y0+h*k3); result=y0+h*(k1+2*k2+2*k3+k4)/6; else double x1,y1; x1=xn-h; y1=Runge_Kuta(f, x0, y0, xn-h,step-1); k1=f(x1,y1); k2=f(x1+h/2, y1+h*k1/2); k3=f(x1+h/2, y1+h*k2/2); k4=f(x1+h, y1+h*k3); result=y1+h*(k1+2*k2+2*

4、k3+k4)/6; return(result);int main()double f(double x, double y); double x0,y0; double a,b;/ int step; cout<<"请输入初值x0,y0:" cin>>x0>>y0; cout<<"请输入区间：" cin>>a>>b; /double x0=0,y0=1; double x,y,step; int i; cout<<"请输入步长：" cin>&g

5、t;step; /step=0.1; cout.precision(10); for(i=0;i<=(b-a)/step;i+)x=x0+i*step; cout<<setw(8)<<x<<setw(18)<<Runge_Kuta(f,x0,y0,x,i)<<endl; return(0); double f(double x, double y) double r; r=(x-y)/2; return(r); 2、 高斯列主元法解线性方程组2.1、 算法说明首先将线性方程组做成增光矩阵，对增广矩阵进行行变换。对第元素，在第i列

6、中，第i行及以下的元素选取绝对值最大的元素，将该元素最大的行与第i行交换，然后采用高斯消元法将新得到的消去第i行以下的元素。一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作，即可以得到方程组的解。2.2、 高斯列主元算法流程图i<=0;j<=0NY开始输入未知数个数mi<m*(m+1)max j<=i:max<=|*(head+(m+1)*i+i)|;k<=i+1|*(head+(m+1)*k+i)|>maxmax <=|*(head+(m+1)*k+i)|;max i<=k;NYYK<mK<=k+1;YNma

7、x i!=itemp<=*(head*(m+1)*i+k);*(head*(m+1)*i+k)<=*(head*(m+1)*max i+k)*(head*(m+1)*max i+k)=tempk<=0Nk<=m+1Y NYi!=jk<=0*(head*(m+1)*i+k)<=*(head+(m+1)*j+k)*(head*(m+1)*j+i) *(head*(m+1)* i+k)/(*head+(m+1)*i+i)k<=m+1YNj<=j+1Yj<mYNi=0NY结束*(head+(m+1)*i+m)<=*(head+(m+1)*i+

8、m)/(*head+(m+1)*i+m)；i=i+1；图2-1 算法流程图输出*(head+(m+1)*i+m)j<m2.3、 高斯列主元程序调试对所编写的高斯列主元程序进行编译和链接，然后执行得如下所示的窗口，我们按命令输入增广矩阵的行数为3，输入3行5列的增广矩阵,运行界面为：图2-2 高斯列主元程序调试2.4、 高斯列主元算法代码#include<iostream> #include<cmath>using namespace std;void load();const N=20;float aNN;int m;int main()int i,j;int c

9、,k,n,p,r;float xN,lNN,s,d;cout<<"下面请输入未知数的个数m="cin>>m;cout<<endl;cout<<"请按顺序输入增广矩阵a:"<<endl;load();for(i=0;i<m;i+) for(j=i;j<m;j+) c=(fabs(aji)>fabs(aii)?j:i; /*找列最大元素*/for(n=0;n<m+1;n+) s=ain; ain=acn; acn=s; /*将列最大数防在对角线上*/for(p=0;p<

10、m+1;p+)cout<<aip<<"t"cout<<endl;for(k=i+1;k<m;k+) lki=aki/aii; for(r=i;r<m+1;r+) /*化成三角阵*/ akr=akr-lki*air; xm-1=am-1m/am-1m-1;for(i=m-2;i>=0;i-) d=0;for(j=i+1;j<m;j+)d=d+aij*xj;xi=(aim-d)/aii;cout<<"该方程组的解为:"<<endl;for(i=0;i<m;i+)cout

11、<<"x"<<i<<"="<<xi<<"t" /system("pause"); return 0; void load()int i,j;for(i=0;i<m;i+)for(j=0;j<m+1;j+)cin>>aij;3、 牛顿法解非线性方程组3.1、 算法说明设已知。第1步：计算函数 （3-1）第2步：计算雅可比矩阵 （3-2）第3步：求线性方程组 的解。第4步：计算下一点 （3-3）重复上述过程。3.2、 牛顿法解非线性方程

12、组算法流程图图3-1 算法流程图3.3、 牛顿法解非线性方程组算法程序调试图3-2牛顿法解非线性方程组算法程序调试应用本程序解方程组， 初始近似值x0,y0分别为2.00和0.25，经过3次迭代求出X(1)=1.900691和X(2)=0.311213。图3-2牛顿法解非线性方程组算法程序运行结果3.4、 牛顿法解非线性方程组算法程序代码#include<iostream>#include<cmath>#define N 2 / 非线性方程组中方程个数、未知量个数 #define epsilon 0.0001 / 差向量1范数的上限#define max 10 /最大迭

13、代次数using namespace std;const int N2=2*N;int main()void ff(float xxN,float yyN);void ffjacobian(float xxN,float yyNN);void inv_jacobian(float yyNN,float invNN);void newdim(float x0N, float invNN,float y0N,float x1N);float x0N=2.0,0.25,y0N,jacobianNN,invjacobianNN,x1N,errornorm;int i,iter=0;cout<&l

14、t;"初始解向量："<<endl;for (i=0;i<N;i+) cout<<x0i<<" " cout<<endl;do iter=iter+1;cout<<"第 "<<iter<<" 次迭代开始："<<endl;/jis uan xiang liang han shu zhi-yin bian liang xiang liang y0ff(x0,y0);/jis uan jacobian ju zhen ja

15、cobianffjacobian(x0,jacobian);/jis uan jacobian ju zhen de ni juzhen invjacobianinv_jacobian(jacobian,invjacobian);/you jie xiang liang x0 ji suan jie xiang liang x1newdim(x0, invjacobian,y0,x1);/ji suan cha xiang liang de 1 fan shuerrornorm=0;for (i=0;i<N;i+) errornorm=errornorm+fabs(x1i-x0i);if

16、 (errornorm<epsilon) break;for (i=0;i<N;i+)x0i=x1i; while (iter<max);return 0;void ff(float xxN,float yyN) float x,y; int i; x=xx0;y=xx1;/非线性方程组yy0=x*x-2*x-y+0.5;yy1=x*x+4*y*y-4; cout<<"因变量向量："<<endl;for( i=0;i<N;i+) cout<<yyi<<" " cout<<

17、endl; cout<<endl;void ffjacobian(float xxN,float yyNN)float x,y;int i,j;x=xx0;y=xx1;yy00=2*x-2;yy01=-1;yy10=2*x;yy11=8*y;cout<<"雅克比矩阵： "<<endl;for( i=0;i<N;i+)for(j=0;j<N;j+) cout<<yyij<<" " cout<<endl;cout<<endl;void inv_jacobian(f

18、loat yyNN,float invNN)float augNN2,L; int i,j,k; cout<<"计算雅克比矩阵的逆: "<<endl; for (i=0;i<N;i+)for(j=0;j<N;j+) augij=yyij; for(j=N;j<N2;j+) if(j=i+N) augij=1; else augij=0; for (i=0;i<N;i+)for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl;cout

19、<<endl;for (i=0;i<N;i+) for (k=i+1;k<N;k+)L=-augki/augii; for(j=i;j<N2;j+) augkj=augkj+L*augij;for (i=0;i<N;i+)for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl;cout<<endl; for (i=N-1;i>0;i-) for (k=i-1;k>=0;k-)L=-augki /augii; for(j=N2-1;j>

20、=0;j-) augkj=augkj+L*augij;for (i=0;i<N;i+)for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl;cout<<endl;for (i=N-1;i>=0;i-)for(j=N2-1;j>=0;j-)augij=augij/augii;for (i=0;i<N;i+)for(j=0;j<N2;j+) cout<<augij<<" " cout<<endl; for

21、(j=N;j<N2;j+) invij-N=augij;cout<<endl;cout<<"雅克比矩阵的逆： "<<endl;for (i=0;i<N;i+) for(j=0;j<N;j+) cout<<invij<<" " cout<<endl; cout<<endl;void newdim(float x0N, float invNN,float y0N,float x1N)int i,j;/计算非线性方程组的近似解向量float sum=0;for

22、(i=0;i<N;i+) sum=0; for(j=0;j<N;j+) sum=sum+invij*y0j;x1i=x0i-sum; cout<<"近似解向量："<<endl;for (i=0;i<N;i+) cout<<x1i<<" "<<endl; 4、 龙贝格求积分算法4.1、 算法说明生成的逼近表，并以为最终解来逼近积分 （4-1）逼近存在于一个特别的下三角矩阵中，第0列元素用基于个a,b子区间的连续梯形方法计算，然后利用龙贝格公式计算。当时，第行的元素为 （4-2）当

23、时，程序在第行结束。4.2、 龙贝格求积分算法流程图图4-1 算法流程图图4-1 算法流程图4.3、 龙贝格求积分算法程序调试我们以求解积分，精度为0.0001，最高迭代10次为例，对所编写的龙贝格求积分算法程序进行编译和链接，经执行后得如下所示的窗口图4-2龙贝格求积分算法程序调试说明：应用Romberg算法求在区间上的精度为0.0001的积分为0.494508。4.4、 龙贝格求积分算法代码#include<iostream>#include<cmath>using namespace std;#define f(x) sin(x*x) /举例函数#define e

24、psilon 0.0001 /精度#define MAXREPT 10 /迭代次数,到最后仍达不到精度要求,则输出T(m=10).double Romberg(double aa, double bb) /aa,bb 积分上下限 int m, n;/m控制迭代次数, 而n控制复化梯形积分的分点数. n=2m double h, x; double s, q; double ep; /精度要求 double *y = new doubleMAXREPT;/为节省空间,只需一维数组 /每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新double p;/p总是指示待计算元素

25、的前一个元素(同一行) /迭代初值 h = bb - aa; y0 = h*(f(aa) + f(bb)/2.0; m = 1; n = 1; ep = epsilon + 1.0; /迭代计算 while (ep >= epsilon) && (m < MAXREPT) /复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增 p = 0.0; for (int i=0; i<n; i+)/求Hn x = aa + (i+0.5)*h; p = p + f(x); p = (y0 + h*p)/2.0;/求T2n = 1/2(Tn+Hn

26、),用p指示 /求第m行元素,根据Romberg计算表本行的前一个元素(p指示), /和上一行左上角元素(yk-1指示)求得. s = 1.0; for (int k=1; k<=m; k+) s = 4.0*s; q = (s*p - yk-1)/(s - 1.0); yk-1 = p; p = q; p = fabs(q - ym-1); m = m + 1; ym-1 = q; n = n + n; h = h/2.0; return (q);int main() double a,b; cout<<"Romberg积分,请输入积分范围a,b:"&l

27、t;<endl; cin>>a>>b; cout<<"积分结果:"<<Romberg(a, b)<<endl; return 0;5、 三次样条插值算法5.1 三次样条插值算法说明表表5-1 三次样条插值算法说明表策略描述包含和的方程(i)三次紧压样条，确定，（如果导数已知，这是“最佳选择”）(ii)natural三次样条（一条“松弛曲线”），(iii)外挂到端点(iv) 是靠近端点的常量，(v)在每个端点处指定，5.2、 三次样条插值算法（压紧样条）程序调试我们将所编写的程序三次样条插值算法（压紧样条）程序

28、进行调试图5-1三次样条插值算法（压紧样条）程序输入界面、运行结果图5-2三次样条插值算法程序运行结果（a） 图5-2三次样条插值算法程序运行结果 （b）运行结果分析： （5-1）作图程序（Matlab）：x1=0:0.01:1;y1=0.06*(x1 - 1).3 + 0.42*(x1 - 0).3 - 0.06*(x1 - 1) + 0.08*(x1 - 0);x2=1:0.01:2;y2=-0.42*(x2 - 2).3 - 0.62*(x2 - 1).3 - 0.08*(x2 - 2) + 2.62*(x2 - 1);x3=2:0.01:3;y3=0.62*(x3 - 3).3 + 0

29、.06*(x3 - 2).3 - 2.62*(x3 - 3) + 1.44*(x3 - 2);X=0 1 2 3;Y=0 0.5 2 1.5;plot(x1,y1,x2,y2,x3,y3,X,Y,'*')gtext('S1')gtext('S2')gtext('S3')图形为：图5-3 三次样条插值算法Matlab作图分析5.3、 三次样条插值算法（压紧样条）代码#include<iostream>#include<iomanip>using namespace std;const int max = 50

30、;float xmax, ymax, hmax;float cmax, amax, fxymmax;float f(int x1, int x2, int x3) float a = (yx3 - yx2) / (xx3 - xx2); float b = (yx2 - yx1) / (xx2 - xx1); return (a - b)/(xx3 - xx1); /求差分void cal_m(int n) /用追赶法求解出弯矩向量M float Bmax; B0 = c0 / 2; for(int i = 1; i < n; i+) Bi = ci / (2 - ai*Bi-1); f

31、xym0 = fxym0 / 2; for(i = 1; i <= n; i+) fxymi = (fxymi - ai*fxymi-1) / (2 - ai*Bi-1); for(i = n-1; i >= 0; i-) fxymi = fxymi - Bi*fxymi+1;void printout(int n);int main() int n,i; char ch; do cout<<"输入x的最大下标:" cin>>n; for(i = 0; i <= n; i+) cout<<"Please put

32、 in X"<<i<<':' cin>>xi; cout<<"Please put in Y"<<i<<':' cin>>yi; for(i = 0; i < n; i+) /求步长 hi = xi+1 - xi; cout<<"Please 输入边界条件n 1 :已知两端的一阶导数n 2 :两端的二阶导数已知n 默认:自然边界条件n" int t; float f0, f1; cin>>t; swi

33、tch(t) case 1:cout<<"输入 Y0' Y"<<n<<"'n" cin>>f0>>f1; c0 = 1; an = 1; fxym0 = 6*(y1 - y0) / (x1 - x0) - f0) / h0; fxymn = 6*(f1 - (yn - yn-1) / (xn - xn-1) / hn-1; break; case 2:cout<<"输入 Y0" Y"<<n<<""

34、;n" cin>>f0>>f1; c0 = an = 0; fxym0 = 2*f0; fxymn = 2*f1; break; default:cout<<"不可用n"/待定 ; for(i = 1; i < n; i+) fxymi = 6 * f(i-1, i, i+1); for(i = 1; i < n; i+) ai = hi-1 / (hi + hi-1); ci = 1 - ai; an = hn-1 / (hn-1 + hn); cal_m(n); cout<<"n输出三次样条

35、插值函数：n" printout(n); cout<<"Do you have another 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<<": "<<xi<<" , "<

36、;<xi+1<<"n"<<"t"cout<<"S"<<i+1<<"=" float t = fxymi/(6*hi); if(t > 0) cout<<-t<<"*(x - "<<xi+1<<")3" else cout<<-t<<"*(x - "<<xi+1<<")3"

37、 t = fxymi+1/(6*hi); if(t > 0) cout<<" + "<<t<<"*(x - "<<xi<<")3" else cout<<" - "<<t<<"*(x - "<<xi<<")3" cout<<"nt" t = (yi - fxymi*hi*hi/6)/hi; if(t > 0) co

38、ut<<"- "<<t<<"*(x - "<<xi+1<<")" else cout<<"- "<<-t<<"*(x - "<<xi+1<<")" t = (yi+1 - fxymi+1*hi*hi/6)/hi; if(t > 0) cout<<" + "<<t<<"*(x - &quo

39、t;<<xi<<")" else cout<<" - "<<-t<<"*(x - "<<xi<<")" cout<<endl; cout<<endl;6、M次多项式曲线拟合6.1、算法说明设有N个点，横坐标是确定的。最小二乘抛物线的系数表示为 （6-1）求解A，B和C的线性方程组为 （6-2）6.2、 M次多项式曲线拟合算法流程图图6-1 算法流程图6.3、 M次多项式曲线拟合算法程序调试我们按命令依次输入命

40、令如下命令后，得程序执行结果如下图6-2 算法调试图作图程序： 图形为：X=-3 -1 1 3;Y=15 5 1 5;x=-3:0.01:3;y=0.875*x.2-1.7*x+2.145;plot(X,Y,'go',x,y,'r-')gtext('拟合曲线')图6-3 图形分析6.4、 M次多项式曲线拟合算法代码#include<iostream>#include<cmath>#define MAX 20using namespace std;/求解任意可逆矩阵的逆，X为待求解矩阵，E为全零矩阵，非单位矩阵，也可以是单位

41、矩阵void inv(double XMAXMAX,int n,double EMAXMAX)int i,j,k;double temp=0;for(i=0;i<MAX;i+)for(j=0;j<MAX;j+)if(i=j)Eij=1;for(i=0;i<n-1;i+)temp=Xii;for(j=0;j<n;j+)Xij=Xij/temp;Eij=Eij/temp;for(k=0;k<n;k+)if(k=i)continue;temp=-Xii*Xki;for(j=0;j<n;j+)Xkj=Xkj+temp*Xij;Ekj=Ekj+temp*Eij;int

42、 main()int n,M,i,j,k;double XMAX=0,YMAX=0,FMAXMAX=0,BMAX=0;double AMAXMAX=0,BFMAXMAX=0,CMAX=0,EMAXMAX=0;cout<<"M次多项式曲线拟合n请先输入待拟合的点的个数:"cin>>n;cout<<"n请输入"<<n<<"个点的X坐标序列:n"for(i=0;i<n;i+)cin>>Xi;cout<<"n请输入"<<n

43、<<"个点的Y坐标序列:n"for(i=0;i<n;i+)cin>>Yi;cout<<"n请输入需要拟合的次数:"cin>>M;for(i=0;i<n;i+)for(k=1;k<=M+1;k+)Fik-1=pow(Xi,k-1);/求F的转置for(i=0;i<n;i+)for(j=0;j<M+1;j+)BFji=Fij; /计算其转置的BF与F的乘 for(i=0;i<M+1;i+)for(j=0;j<M+1;j+)for(k=0;k<n;k+)Aij+=B

44、Fik*Fkj;/计算F的转置BF与Y的乘for(i=0;i<M+1;i+)for(j=0;j<n;j+)Bi+=BFij*Yj;/调用inv函数求解矩阵A的逆矩阵Einv(A,n,E);/计算A的逆BF与B的乘for(i=0;i<M+1;i+)for(j=0;j<n;j+)Ci+=Eij*Bj;cout<<"n拟合后的"<<M<<"次多项式系数为，幂次由高到低：n"for(i=M;i>=0;i-)cout<<Ci<<"t"cout<<

45、;"n拟合后的"<<M<<"次多项式为:n"cout<<"nP(x)="for(i=M;i>=0;i-)if(i=0)cout<<"+"<<Ci;elsecout<<"+"<<Ci<<"*x"<<i;cout<<endl;return 0;7、 牛顿-拉弗森迭代解非线性方程7.1、 牛顿-拉弗森迭代解非线性方程算法概要使用初始近似值，利用迭代 其中计算

46、函数的根的近似值。7.2、 牛顿-拉弗森迭代解非线性方程算法程序调试图7-1牛顿-拉弗森迭代解非线性方程算法程序调试7.3、 牛顿-拉弗森迭代解非线性方程算法代码#include<iostream>using namespace std;#include<cmath>double f(double x)return x*x-2*x-1double df(double x)return 2*x-2;int main() int k=0,max1=0; double p0,delta,epsilon,p1,err,relerr,y;cout<<"牛顿拉

47、弗森法解非线性方程f(x)=x2-2x-1"<<endl;cout<<"初始值为p0=4.8"<<endl; p0=4.8,delta=0.0001,epsilon=0.0001,max1=100;for(k=1;k<=max1;k+)p1=p0-f(p0)/df(p0);err=fabs(p1-p0);relerr=2*err/(fabs(p1)+delta);p0=p1;y=f(p0);if(err<delta)|(relerr<delta)|(fabs(y)<epsilon)break;cout&l

48、t;<"方程的根的近似值为:"<<p0<<endl;cout<<"方程的根的误差估计为:"<<err<<endl;cout<<"迭代次数为:"<<k<<endl;cout<<"方程在p0点的函数值f(p0):"<<y<<endl;return 0; 8、 不动点法解非线性方程8.1、 算法说明先将改写成然后对进行迭代，即 其中然后判断是否成立，成立则返回，不成立就重复以上步骤8.

49、2、 不动点法解非线性方程算法程序调试我们将编写好的不动点迭代法解非线性方程算法程序进行编译，链接和执行后得如下所示结果。图8-1 程序调试图8.3、 不动点法解非线性方程算法代码#include<iostream>#include<cmath>#include<iomanip>#define MAX 20#define eps 1e-10using namespace std;double g(double x)return 1-x*x/4+x; main()double PMAX=0,err=0.0,relerr=0.0,tol=0.0,p=0.0,p0

50、=0.0;int k=0,max1=0,i=0;cout<<"不动点法解非线性方程f(x)=1-x2/2"<<endl;cout<<"方程在0，1上有解，初始值为p0=0"<<endl;/初始化P0=p0=0;max1=100;tol=0.001;for(k=2;k<=max1;k+)Pk-1=g(Pk-2);err=fabs(Pk-1-Pk-2);relerr=err/(fabs(Pk-1+eps);p=Pk-1;if(err<tol)|(relerr<tol)break;if(k=ma

51、x1)cout<<"迭代次数超过允许的最大迭代次数！"<<endl;cout<<"不动点的近似值为:"<<p<<endl; cout<<"程序迭代次数为:"<<k<<endl;cout<<"近似值的误差为:"<<err<<endl;cout<<"求解不动点近似值的序列:"<<endl;for(i=0;i<k;i+)cout<<

52、;setw(10)<<Pi;cout<<endl;return 0;9. 拉格朗日插值9.1、 算法说明对每个固定的k,拉格朗日系数多项式具有性质为：. （9-1）9.2、 拉格朗日插值算法程序调试首先编写好程序，然后编译链接，运行程序，按照程序提示依次输入点的个数、点数对应的x值、y值、纵坐标、再输入x值，输入结果如下：图9-1拉格朗日插值算法程序调试图9-2拉格朗日插值算法程序计算结果9.3、 拉格朗日插值算法程序代码：#include<iostream>using namespace std;double func(double X,int k,dou

53、ble x,int n);int main()double Sn=0; int n;cout<<"请输入点的个数n:"cin>>n;double*x=(double*)malloc(n*sizeof(double); double*y=(double*)malloc(n*sizeof(double);double X;int i;for(i=0;i<n;i+)cout<<"请输入x"<<i+1<<",y"<<i+1<<":"&

54、lt;<endl;cin>>xi>>yi;cout<<"请输入x"cin>>X;for(i=0;i<n;i+) Sn=Sn+func(X,i,x,n)*yi;cout<<"通过拉格朗日插值公式所得的结果："<<"当x="<<X<<"时，y="<<Sn<<endl; return 0;double func(double X,int k,double x,int n)int i;double Pn=1;double p;for(i=0;i<n;i+)if(i=k) continue