计算方法实验报告.docx
- 文档编号:6733665
- 上传时间:2023-01-09
- 格式:DOCX
- 页数:22
- 大小:105.39KB
计算方法实验报告.docx
《计算方法实验报告.docx》由会员分享,可在线阅读,更多相关《计算方法实验报告.docx(22页珍藏版)》请在冰豆网上搜索。
计算方法实验报告
数值分析与算法
实
验
报
告
学院:
指导老师:
班级:
姓名:
学号:
实习题一
1.用两种不同的顺序计算
,分析其误差的变化。
●实验代码
#include
#include
usingnamespacestd;
voidmain(){
floatnum=0;
inti;
for(i=0;i<10000;i++){
num=num+1.0/((i+1)*(i+1));
}
cout<<"从小加到大的结果:
"< num=0; for(i=9999;i>=0;i--){ num=num+1.0/((i+1)*(i+1)); } cout<<"从大加到小的结果: "< } ●运行窗口 4.设 ,已知其精确值为1/2(3/2-1/N-1/N+1)。 1)编制按从大到小的顺序计算Sn的程序; 2)编制按从小到大的顺序计算Sn的程序; 3)按两种顺序分别计算S1000,S10000,S30000,并指出有效位数。 ●实验代码 #include #include usingnamespacestd; doublenum; doublesuml(intN){ num=0; intj; for(j=2;j num=num+1.0/(j*j-1); } returnnum; } doublesumh(intN){ num=0; intj; for(j=N;j>1;j--){ num=num+1.0/(j*j-1); } returnnum; } voidmain(){ intN; cout<<"inputN"< cin>>N; cout<<"从小到大加的结果: "< cout<<"从大到小加的结果: "< } ●运行窗口 实验1 拉格朗日插值法 一、方法原理 n次拉格朗日插值多项式为: Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x) n=1时,称为线性插值,L1(x)=y0(x-x1)/(x0-x1)+y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0) n=2时,称为二次插值或抛物线插值,精度相对高些 L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1) 二、主要思路 使用线性方程组求系数构造插值公式相对复杂,可改用构造方法来插值。 对节点xi(i=0,1,…,n)中任一点xk(0<=k<=n)作一n次多项式lk(xk),使它在该点上取值为1,而在其余点xi(i=0,1,…,k-1,k+1,…,n)上为0,则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+…+ynln(x) 上式表明: n个点xi(i=0,1,…,k-1,k+1,…,n)都是lk(x)的零点。 可求得lk 三.计算方法及过程: 1.输入节点的个数n 2.输入各个节点的横纵坐标 3.输入插值点 4.调用函数,返回z 函数语句与形参说明 程序源代码如下: 形参与函数类型 参数意义 intn 节点的个数 doublex[n](double*x) 存放n个节点的值 doubley[n](double*y) 存放n个节点相对应的函数值 doublep 指定插值点的值 doublefun() 函数返回一个双精度实型函数值,即插值点p处的近似函数值 #include #include usingnamespacestd; #defineN100 doublefun(double*x,double*y,intn,doublep); voidmain() {inti,n; cout<<"输入节点的个数n: "; cin>>n; doublex[N],y[N],p; cout<<"pleaseinputxiangliangx="< for(i=0;i cout<<"pleaseinputxiangliangy="< for(i=0;i cout<<"pleaseinputLagelangrichazhiJieDianp="< cin>>p; cout<<"TheAnswer="< system("pause");} doublefun(doublex[],doubley[],intn,doublep) {doublez=0,s=1.0; intk=0,i=0; doubleL[N]; while(k {if(k==0) {for(i=1;i L[0]=s*y[0]; k=k+1;} else {s=1.0; for(i=0;i<=k-1;i++)s=s*((p-x[i])/(x[k]-x[i])); for(i=k+1;i L[k]=s*y[k]; k++;} } for(i=0;i returnz; } 四.运行结果测试: 五.实验分析 n=2时,为一次插值,即线性插值 n=3时,为二次插值,即抛物线插值 n=1,此时只有一个节点,插值点的值就是该节点的函数值 n<1时,结果都是返回0的;这里做了n=0和n=-7两种情况 3 常用的是线性插值和抛物线插值,显然,抛物线精度相对高些 n次插值多项式Ln(x)通常是次数为n的多项式,特殊情况可能次数小于n.例如: 通过三点的二次插值多项式L2(x),如果三点共线,则y=L2(x)就是一条直线,而不是抛物线,这时L2(x)是一次式。 拟合曲线光顺性差 实验2 牛顿插值法 一、方法原理及基本思路 在拉格朗日插值方法中,若增加一个节点数据,其插值的多项式需重新计算。 现构造一个插值多项式Nn(x),只需对Nn-1(x)作简单修正(如增加某项)即可得到,这样计算方便。 利用牛顿插值公式,当增加一个节点时,只需在后面多计算一项,而前面的计算仍有用;另一方面Nn(x)的各项系数恰好又是各阶差商,而各阶差商可用差商公式来计算。 由线性代数知,对任何一个不高n次的多项式P(x)=b0+b1x+b2x2+…+bnxn(幂基)① 也可将其写成P(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)+…+an(x-x0)…(x-xn-1) 其中ai为系数,xi为给定节点,可由①求出ai 一般情况下,牛顿插值多项式Nn(x)可写成: Nn(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)+…+an(x-x0)…(x-xn-1)) 只需求出系数ai,即可得到插值多项式。 二、计算方法及过程 1.先后输入节点个数n和节点的横纵坐标,插值点的横坐标,最后输入精度e 2.用do-while循环语句得到跳出循环时k的值 3.将k值与n-1进行比较,若在达到精度时k 函数语句与形参说明 形参与函数类型 参数意义 intn 节点的个数 floatx[MAX] 存放n个节点的值(从小到大) Floaty[MAX]; 存放n个节点相对应的函数值 floatx0,y0 指定插值点的横纵坐标 floate 精度 程序源代码如下: #include #include usingnamespacestd; #defineMAX100 voidmain() { floatx[MAX],y[MAX]; floatx0,y0,e,N1; floatN0=0; inti,k,n; cout<<"输入节点的个数,n="; cin>>n; cout<<"请输入插值节点! "< for(i=0;i cout<<"请输入对应的函数值! "< for(i=0;i cout<<"请输入插值点,x0="; cin>>x0; cout<<"请输入精度,e="; cin>>e; y0=1; N1=y[0]; k=0; do { k++; N0=N1; y0=y0*(x0-x[k-1]); for(i=0;i N1=N0+y0*y[k]; }while(fabs(N1-N0)>e&&k<(n-1)); if(k==(n-1))cout<<"计算失败! "; if(k<(n-1))cout<<"输出结果y[x0]="< system("pause"); } 三.运行结果测试: 题一: 已知f(x)=sh(x)的函数表如下: 计算f(0.23)的近似值 xi 0 0.20 0.30 0.50 Sh(xi) 0 0.20134 0.30452 0.52110 题二: 已知根号100等于10,根号121等于11,根号144等于12; 用二次插值和三次插值计算根号115的近似值。 四.实验分析 显然,题二中当n=2,n=3时,k的值没有在n-1前满足各自的精度要求,所以计算失败 从题一中可以看出,插值点的个数、精度、插值点的选择都会影响实验的结果; 我们通常会选择与插值点最接近的节点,可以提高精度; 在可以计算出结果的情况下,插值点越多,结果越精确 实习题二 1.用牛顿法求下列方程的根: 1) ●实验代码 #include #include #include usingnamespacestd; doublefy(doublex){ returnx*x-exp(x); } doublefd(doublex){ return2*x-exp(x); } voidmain(){ inti=0; doubley=2.0; for(i=0;i<50;i++){ y=y-fy(y)/fd(y); cout< i++; } } ●运行窗口 2) 原理同1 (1),故略 2.编写一个割线法程序,求解上述个方程 编写1 (1)的割线法程序 ●实验代码 #include #include usingnamespacestd; doublef(doublex){ returnx*x-exp(x); } doublefd(doublex){ return2*x-exp(x); } voidmain(){ inti=0; doublex1=0,x2=1,x3;//该方程在大于0上没根,exp(x)的增长速度比x*x大的多 while (1){ x3=x2-(f(x2)/(f(x2)-f(x1)))*(x2-x1); x1=x2; x2=x3; cout< i++; if(i>50)break; } } ●运行窗口 实习题三 1.用列主元消去法解方程组: 1) ●实验代码 #include #include #include usingnamespacestd; voidmain(){ voidColPivot(float*,int,float[]); inti; floatx[3]; floatc[4][5]={1,1,0,3,4,2,1,-1,1,1,3,-1,-1,3,-3,-1,2,3,-1,4}; ColPivot(c[0],4,x); for(i=0;i<=3;i++)printf("x[%d]=%f\n",i,x[i]); } voidColPivot(float*c,intn,floatx[]){ inti,j,k,t; floatp; for(i=0;i<=n-2;i++){ k=i; for(j=i+1;j if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))k=j; if(k! =i) for(j=i;j<=n;j++){ p=*(c+i*(n+1)+j); *(c+i*(n+1)+j)=*(c+k*(n+1)+j); *(c+k*(n+1)+j)=p; } for(j=i+1;j p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i)); for(t=i;t<=n;t++)*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t)); } } for(i=n-1;i>=0;i--){ for(j=n-1;j>=i+1;j--) (*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j)); x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i)); } } ●运行窗口 4.编写用追赶法解三对角线性方程组的程序,并解下列方程组: 2) ●实验代码 #include #include usingnamespacestd; voidmain(){ floata=1; floatb=-4; floatc=1; floatd[11]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15}; floatl[11]; floatbb[11]; floaty[11]; floatx[11]; bb[1]=b; y[1]=d[1]; inti; for(i=2;i<11;i++){ l[i]=a/bb[i-1]; bb[i]=b-l[i]*c; y[i]=d[i]-l[i]*y[i]; } x[10]=y[10]/bb[10]; for(i=9;1>0;i--){ x[i]=(y[i]-c*x[i+1])/bb[i]; } for(i=1;i<11;i++) cout<<'x'< '< } ●运行窗口 #include #include #include #include #include #defineN11 main() { floata[N]={0,0,1,1,1,1,1,1,1,1,1}; floatb[N]={0,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4}; floatc[N]={0,1,1,1,1,1,1,1,1,1,0}; floatd[N]={0,-27,-15,-15,-15,-15,-15,-15,-15,-15,-15}; floatx[N]={0,0,0,0,0,0,0,0,0,0,0}; floatr[N]={0,0,0,0,0,0,0,0,0,0,0}; floaty[N]={0,0,0,0,0,0,0,0,0,0,0}; floatq; //clrscr(); intk; r[1]=c[1]/b[1]; y[1]=d[1]/b[1]; for(k=2;k { q=b[k]-r[k-1]*a[k]; r[k]=c[k]/q; y[k]=(d[k]-y[k-1]*a[k])/q; } y[N-1]=(d[N-1]-y[N-2]*a[N-1])/(b[N-1]-r[N-2]*a[N-1]); x[N-1]=y[N-1]; for(k=N-2;k>=1;k--) x[k]=y[k]-r[k]*x[k+1]; for(k=1;k printf("x[%d]=%f\n",k,x[k]); getch(); return0; } ●运行窗口 5. 雅克比迭代法 ●实验代码 #include #include #defineeps1e-6 #definemax100 voidJacobi(float*a,intn,floatx[]) { inti,j,k=0; floatepsilon,s; float*y=newfloat[n]; for(i=0;i while (1) { epsilon=0; k++; for(i=0;i { s=0; for(j=0;j { if(j==i)continue; s+=*(a+i*(n+1)+j)*x[j]; } y[i]=(*(a+i*(n+1)+n)-s)/(*(a+i*(n+1)+i)); epsilon+=fabs(y[i]-x[i]); } for(i=0;i if(epsilon {printf("迭代次数为: %d\n",k);return;} if(k>=max) {printf("迭代发散");return;} } deletey; } voidmain() { inti; floata[4][5]={10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25}; floatx[4]; Jacobi(a[0],4,x); for(i=0;i<4;i++)printf("x[%d]=%f\n",i,x[i]); } ●运行窗口 高斯-塞德尔迭代法 ●实验代码 #include #include #defineN500 voidmain() { inti; floatx[4]; floatc[4][5]={10,-1,2,0,-11,0,8,-1,3,-11,2,-1,10,0,6,-1,3,-1,11,25}; voidGaussSeidel(float*,int,float[]); GaussSeidel(c[0],4,x); for(i=0;i<4;i++)printf("x[%d]=%f\n",i,x[i]); } voidGaussSeidel(float*a,intn,floatx[]) { inti,j,k=1; floatd,dx,eps; for(i=0;i while (1) { eps=0; for(i=0;i { d=0; for(j=0;j { if(j==i)continue; d+=*(a+i*(n+1)+j)*x[j]; } dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i)); eps+=fabs(dx-x[i]); x[i]=dx; } if(eps<1e-6){printf("迭代次数为: %d\n",k);return;} if(k>N) { printf("迭代发散\n"); return; } k++; } } ●运行窗口 实习题四 2.按下列数据 Xi 0.30 0.42 0.50 0.58 0.66 0.72 Yi 1.04403 1.08462 1.11803 1.15603 1.19817 1.23223 作5次插值,并求X1=0.46,X2=0.55,X3=0.60时的函数近似值。 ●实验代码 #include #defineN5 voidDifference(floatx[],floaty[],intn) { float*f=newfloat[n+1]; intk,i; for(k=1;k<=n;k++) { f[0]=y[k]; for(i=0;i f[i+1]=(f[i]-y[i])/(x[k]-x[i]); y[k]=f[k]; } deletef; return; } voidmain() { inti; floata,b,c,varx=0.46,vary=0.55,varz=0.60; floatx[N+1]={0.30,0.42,0.50,0.58,0.66,0.72}; floaty[N+1]={1.04403,1.08462,1.11803,1.15603,1.19817,1.23223}; Difference(x,y,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算方法 实验 报告
![提示](https://static.bdocx.com/images/bang_tan.gif)