计算方法上机作业
西安交通大学计算方法B大作业
计算方法上机报告姓名：学号：班级：目录题目一------------------------------------------------------------------------------------------ - 4 -1.1题目内容 ---------------------------------------------------------------------------- - 4 -1.2算法思想 ---------------------------------------------------------------------------- - 4 -1.3Matlab源程序----------------------------------------------------------------------- - 5 -1.4计算结果及总结 ------------------------------------------------------------------- - 5 - 题目二------------------------------------------------------------------------------------------ - 7 -2.1题目内容 ---------------------------------------------------------------------------- - 7 -2.2算法思想 ---------------------------------------------------------------------------- - 7 -2.3 Matlab源程序---------------------------------------------------------------------- - 8 -2.4计算结果及总结 ------------------------------------------------------------------- - 9 - 题目三----------------------------------------------------------------------------------------- - 11 -3.1题目内容 --------------------------------------------------------------------------- - 11 -3.2算法思想 --------------------------------------------------------------------------- - 11 -3.3Matlab源程序---------------------------------------------------------------------- - 13 -3.4计算结果及总结 ------------------------------------------------------------------ - 14 - 题目四----------------------------------------------------------------------------------------- - 15 -4.1题目内容 --------------------------------------------------------------------------- - 15 -4.2算法思想 --------------------------------------------------------------------------- - 15 -4.3Matlab源程序---------------------------------------------------------------------- - 15 -4.4计算结果及总结 ------------------------------------------------------------------ - 16 - 题目五----------------------------------------------------------------------------------------- - 18 -5.1题目内容 --------------------------------------------------------------------------- - 18 -5.2算法思想 --------------------------------------------------------------------------- - 18 -5.3 Matlab源程序--------------------------------------------------------------------- - 18 -5.3.1非压缩带状对角方程组------------------------------------------------- - 18 -5.3.2压缩带状对角方程组---------------------------------------------------- - 20 -5.4实验结果及分析 ------------------------------------------------------------------ - 22 -5.4.1Matlab运行结果 ---------------------------------------------------------- - 22 -5.4.2总结分析------------------------------------------------------------------- - 24 -5.5本专业算例 ------------------------------------------------------------------------ - 24 - 学习感悟-------------------------------------------------------------------------------------- - 27 -题目一1.1题目内容计算以下和式：0142111681848586n n S n n n n ∞=⎛⎫=--- ⎪++++⎝⎭∑，要求： (1)若保留11个有效数字，给出计算结果，并评价计算的算法； (2)若要保留30个有效数字，则又将如何进行计算。
西安电子科技大学出版社计算方法上机答案
西安电子科技大学出版社《计算方法》任传祥等编著第九章计算方法上机参考答案实验一，算法一#include <stdio.h>#include <math.h>double I0=log(6)/log(5),I1;int n=1;main (){while(1){I1=1.0/(n)-I0*5.0;printf("%d %lf\n", n,I1);if(n>=20)break;elseI0=I1;n++;}}实验一，算法二#include <stdio.h>#include <math.h>double I0=(1/105.0+1/126.0)/2,I1;int n=20;main (){printf("%d %lf\n", n,I0);while(1){I1=1.0/(5.0*n)-I0/5.0;printf("%d %lf\n", n-1,I1);if(n<2)break;elseI0=I1;n--;}}实验二，二分法#include <stdio.h>#include <math.h>#define esp 1e-3double f(double x);main (){double a=1,b=2,x;while(fabs(b-a)>esp){x=(a+b)/2;printf("x=%lf\n",x);if(f(x)==0)break;elseif(f(x)*f(a)<0)b=x;elsea=x;}}double f(double x){return pow(x,3)-x-1;}实验二，牛顿迭代法#include<stdio.h>#include<math.h>double f(double x);double f1(double x);#define esp 1e-3void main(){double x0 = 1.5, x1;x1 = x0 - f(x0) / f1(x0);printf("x=%lf\n", x1);x0 = x1;x1 = x0 - f(x0) / f1(x0);printf("x=%lf\n", x1);while (fabs(x1 - x0)>esp){x0 = x1;x1 = x0 - f(x0) / f1(x0);printf("x=%lf\n", x1);} }double f(double x){return pow(x, 3) - x - 1;} double f1(double x){return 3 * x*x - 1;}弦割法#include<stdio.h>#include<math.h>double f(double x);#define esp 1e-3void main(){double x0 = 1.5, x1=2.0,x2;do{ x2=x1 - (x1-x0)*f(x1) /(f(x1)-f(x0));x0=x1;x1=x2;printf("x=%lf\n", x1);}while (fabs(x1 - x0)>esp);{printf("x=%lf\n", x1);}}double f(double x){return pow(x, 3) - x - 1;}实验3#include <stdio.h>/*列主元高斯消去法*/#include <math.h>float x[3],temp,max;float A[3][4]={10,-2,-1,3,-2,10,-1,15,-1,-2,5,10},c[3][4]={10,-2,-1,3,-2,10,-1,15,-1,-2,5,10}; int n=3,i,k,j,m;void main(){for(i=0;i<n;i++){max=A[i][i];k=i;for(j=j+1;j<n;j++){{max=fabs(A[j][i]);k=j;}}if(k!=i){for(j=i+1;j<=n;j++){temp=A[i][j];A[i][j]=A[k][j];A[k][j]=temp;}}for(j=i+1;j<n;j++)for(m=i+1;m<=n;m++){c[j][m]=c[j][m]+(-c[j][i]/c[i][i])*c[i][m];}}for(i=n-1;i>=0;i--){temp=0.0;for(j=n-1;j>=i+1;j--)temp=temp+c[i][j]*x[j];x[i]=(c[i][n]-temp)/c[i][i];}printf("x[1]=%f\nx[2]=%f\nx[3]=%f\n",x[0],x[1],x[2]);实验四，拉格朗日插值#include<stdio.h>int n=5,i,j;double l,L=0,X=0.5;main(){double x[5]={0.4,0.55,0.65,0.8,0.9};doubley[5]={0.41075,0.57815,0.69675,0.88811,1.02652}; for(i=0;i<n;i++){l=y[i];for(j=0;j<n;j++){if(j!=i)l=l*(X-x[j])/(x[i]-x[j]); } L=L+l;}printf("%lf\n",L);return 0;} X=0.5 X=0.7 X=0.85牛顿插值法#include<stdio.h>#include<math.h>main(){double x[5]={0.4,0.55,0.65,0.8,0.9};doubley[5]={0.41075,0.57815,0.69675,0.88811,1.02652};int n=5,i,j;double z;printf("input z\n");scanf("%lf",&z);double a[5][5];for(i=0;i<5;i++)a[i][0]=y[i];for(i=1;i<5;i++)for(j=i;j<5;j++)a[j][i]=(a[j][i-1]-a[j-1][i-1])/(x[j]-x[j-i]);double N=a[0][0],temp=1.0;for(i=1;i<n;i++){temp=temp*(z-x[i-1]);N=N+a[i][i]*temp;}printf("N=%lf\n",N);return 0;}实验五曲线拟合#include <stdio.h>#include <math.h>float x[5]={1,2,3,4,5};float y[5]={7,11,17,27,40};float A[2][3],c[2][3];float z[2],temp,max;int i,j,k,m;int n=2;void main(){for(i=0;i<5;i++){c[0][0]=A[0][0]+=1;c[0][1]=A[0][1]+=x[i];c[0][2]=A[0][2]+=y[i];c[1][0]=A[1][0]+=x[i];c[1][1]=A[1][1]+=x[i]*x[i];c[1][2]=A[1][2]+=x[i]*y[i];}/* for(i=0;i<2;i++){printf(" %lf %lf %lf\n",A[i][0],A[i][1],A[i ][2]);}*/for(i=0;i<n;i++){max=A[i][i];k=i;for(j=j+1;j<n;j++){if(fabs(A[j][i])>max){max=fabs(A[j][i]);k=j;}} if(k!=i){for(j=i+1;j<=n;j++){temp=A[i][j];A[i][j]=A[k][j];A[k][j]=temp;}}for(j=i+1;j<n;j++)for(m=i+1;m<=n;m++){c[j][m]=c[j][m]+(-c[j][i]/c[i][i])*c[i][m];}}for(i=n-1;i>=0;i--){temp=0.0;for(j=n-1;j>=i+1;j--)temp=temp+c[i][j]*z[j];z[i]=(c[i][n]-temp)/c[i][i];}printf("a=%f\nxb=%f\n",z[0],z[1]); }实验六数值积分/*梯形*/#include<stdio.h>#include<math.h> double f(double x); main(){double x[10],y[10];double h,b=1,a=0,I;int n,i;printf("n\n");scanf("%d",&n);h=(b-a)/n;for(i=0;i<=n;i++){x[i]=a+(i*h);y[i]=f(x[i]);}I=f(a)+f(b);for(i=1;i<=n-1;i++){I=I+2*y[i];}I=(h/2)*I;printf("%lf",I);}double f(double x){double f;f=1.0/(1.0+(x*x));return(f);}/*辛普森*/#include<stdio.h>#include<math.h>double f(double x);main(){double x[30],y[30];double h,b=1,a=0,I;int n,i;printf("n\n");scanf("%d",&n);//点乘2扩展h=(b-a)/n;x[10]=1;y[10]=f(x[10]);for(i=0;i<n;i++){x[2*i]=a+(i*h);y[2*i]=f(x[2*i]);x[2*i+1]=a+(i+(1.0/2.0))*h;y[(2*i)+1]=f(x[(2*i)+1]);}I=f(a)+f(b);for(i=0;i<n;i++){I=I+4*y[(2*i)+1];}for(i=1;i<n;i++){I=I+2*y[2*i];}I=(h/6)*I;printf("%lf\n",I);}double f(double x){double f;f=1.0/(1.0+(x*x));return(f);}/*梯形*//*辛普森*/。
东南大学计算方法与实习上机实验一
东南大学计算方法与实习实验报告学院：电子科学与工程学院学号：06A*****姓名：***指导老师：***实习题14、设S N=Σ(1)编制按从大到小的顺序计算S N的程序；(2)编制按从小到大的顺序计算S N的程序；(3)按两种顺序分别计算S1000，S10000，S30000，并指出有效位数。
解析：从大到小时，将S N分解成S N-1=S N-,在计算时根据想要得到的值取合适的最大的值作为首项；同理从小到大时，将S N=S N-1+ ，则取S2=1/3。
则所得式子即为该算法的通项公式。
(1)从大到小算法的C++程序如下：/*从大到小的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;const int max=34000; //根据第(3)问的问题，我选择了最大数为34000作为初值void main(){int num;char jus;double cor,sub;A: cout<<"请输入你想计算的值"<<'\t';cin>>num;double smax=1.0/2.0*(3.0/2.0-1.0/max-1.0/(max+1)),temps;double S[max];// cout<<"s["<<max<<"]="<<setprecision(20)<<smax<<'\n';for(int n=max;n>num;){temps=smax;S[n]=temps;n--;smax=smax-1.0/((n+1)*(n+1)-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/num-1.0/(num+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-smax); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<n<<"]="<<setprecision(20)<<smax<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(2)从小到大算法的C++程序如下：/*从小到大的算法*/#include<iostream>#include<iomanip>#include<cmath>using namespace std;void main(){int max;A: cout<<"请输入你想计算的数，注意不要小于2"<<'\t';cin>>max;double s2=1.0/3.0,temps,cor,sub;char jus;double S[100000];for(int j=2;j<max;){temps=s2;S[j]=temps;j++;s2+=1.0/(j*j-1.0);}cor=1.0/2.0*(3.0/2.0-1.0/j-1.0/(j+1.0)); //利用已知精确值公式计算精确值sub=fabs(cor-s2); //double型取误差的绝对值cout<<"用递推公式算出来的s["<<j<<"]="<<setprecision(20)<<s2<<'\n';cout<<"实际精确值为"<<setprecision(20)<<cor<<'\n';cout<<"则误差为"<<setprecision(20)<<sub<<'\n';cout<<"是否继续计算S[N],是请输入Y，否则输入N！"<<endl;cin>>jus;if ((int)jus==89||(int)jus==121) goto A;}(3)(注：因为程序中setprecision(20)表示输出数值小数位数20，则程序运行时所得到的有效数字在17位左右)ii.选择从小到大的顺序计算S1000、S10000、S30000的值需要计算的项S1000S10000S30000计算值0.74900049950049996 0.74966672220370571 0.74996666722220728实际精确值0.74900049950049952 0.74990000499950005 0.74996666722220373误差 4.4408920985006262*10-16 5.6621374255882984*10-15 3.5527136788005009*10-15有效数字17 17 17附上部分程序运行图：iii.实验分析通过C++程序进行计算验证采用从大到小或者从小到大的递推公式算法得到的数值基本稳定且误差不大。
应用计算方法
结果分析：比较发现，经过两种改进迭代法，求重根时迭代速度明显加快。
3-4 试验目的体验 Steffensen’s method 加速技巧 试验内容：先用 Newton 法求解方程 x-tanx=0 再用 Steffensen 法求解，比较迭代步数。精确到 10-4。 迭代公式： P(k+1)=P(k)-(P(k)-tan(P(k)))/(1-(sec(P（k）))^2) 初值 P0=1 Newton 法： function [ p,k ]=fnewton( p0,max,tol ) for k=1:max p=p0-(p0-tan(p0))/(1-(sec(p0))^2); if abs(p-p0)<tol break; end p0=p; end disp(p); disp(k) % fnewton( 1,100,10^(-4) ) Steffensen 法： function [ p1,k ]=steffensen( p0,max,tol ) n=1; p(1)=p0; while n<=max for k=1:2 p(k+1)=p(k)-(p(k)-tan(p(k)))/(1-(sec(p(k)))^2); end p1=p(1)-(p(2)-p(1))^2/(p(3)-2*p(2)+p(1)); f0=p1-(p1-tan(p1))/(1-(sec(p1))^2); if abs(f0)<tol break; end n=n+1; p(1)=p1; end disp(p1); disp(n) % steffensen( 1,100,10^(-4) ) >> steffensen( 1,100,10^(-4) ) -1.3387e-07 3
3-2 试验目的:考察 Newton 法求单根的收敛速度 应用 Newton 迭代法求 3-1 中的方程，并与 3-1 中的迭代法相比较，考察收敛速度。精确到时 10-4 迭代公式： P(k+1)= P(k)-(2P(k)-eP(k)+3)/(2- eP(k)) 初值 P0=1 和-1。 Newton 法: function [ p,k ] = fnewton( p0,max,tol ) for k=1:max p=p0-(2*p0-exp(p0)+3)/(2-exp(p0)); if abs(p-p0)<tol break; end
计算方法与计算 实验一误差分析
% 输出的量--每次迭代次数k和迭代值xk,
%
--每次迭代的绝对误差juecha和相对误差xiangcha，
误差分析
误差问题是数值分析的基础，又是数值分析中一个困难的课题。在实际计算 中，如果选用了不同的算法，由于舍入误差的影响，将会得到截然不同的结果。 因此，选取算法时注重分析舍入误差的影响，在实际计算中是十分重要的。同时， 由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差，因此算法 的好坏会影响到数值结果的精度。 一、实验目的
因为运行后输出结果为： y 1.370 762 168 154 49， yˆ =1.370 744 664 189
38， R 1.750 396 510 491 47e-005， WU= 1.782 679 830 970 664e-005 104 . 所
以， yˆ 的绝对误差为 10 4 ，故 y
③ 运行后输出计算结果列入表 1–1 和表 1-2 中。
④ 将算法 2 的 MATLAB 调用函数程序的函数分别用 y1=15-2*x^2 和
y1=x-(2*x^2+x-15)/(4*x+1)代替，得到算法 1 和算法 3 的调用函数程序，将其保
存，运行后将三种算法的前 8 个迭代值 x1, x2 ,, x8 列在一起（见表 1-1），进行
的精确解 x* 2.5 比较，观察误差的传播.
算法 1 将已知方程化为同解方程 x 15 2x2 .取初值 x0 2 ，按迭代公式
xk1 15 2xk2
计算方法大作业1 克服Runge现象
x3
x2
x
1
S1 ( x)
-0.34685
0.2086
0.073964
0.038462
S2 (x)
S (xi 0 ) S x(i 0 )
S
'
(xi
0) S
xi' (
0 )i
S
'
'
x(i
0)S
xi' ' (
0)
1 ,n2, . . . , 1
（1）
这里共有了 3n-3 个条件，再加上条件（2）中的 n+1 个插值条件，共有 4n-2 个条件，
因此还需要 2 个方程才能确定 S (x) .通常可在区间[a, b]的端点 a x0,b xn 上各加一个边
dn1
1
2
Mn
dn
（6）
2 1
2
2
2
1 M1 d1
M2
d2
n 1
2
n
1
M
n
1
dn1
n
n 2 M n dn
由式（1）内点拼接条件，可得
i M i1 2M i i M i1 d j i 1, 2,..., n 1
（3） （4）
其中
i
hi 1 hi1
, hi
i
hi hi 1
计算方法上机程序
1.对分+扫描Private Function f(x!)f = x ^ 4 - 5 * x ^ 2 + x + 2 End FunctionPrivate Sub Form_Click() Dim a!, b!, h!, c!, p!, q!, x!a = InputBox("输入a")b = InputBox("输入b")h = InputBox("输入h")x = aDo While x < bIf f(x) * f(x + h) <= 0 Thenp = x: q = x + hDo While Abs(q - p) > 0.00001 c = (p + q) / 2If f(c) = 0 ThenExit DoElseIf f(p) * f(c) < 0 Thenq = cElsep = cEnd IfEnd IfLoopPrint "["; p, q; "]"; cEnd Ifx = x + hLoopEnd Sub2.用牛顿法求a的立方根，精度要求0.000005 Private Sub Form_Click()Dim x0 As Single, x1 As SingleDim a As Integera = InputBox("输入a")If a = 0 ThenPrint "a的立方根=0"EndEnd Ifx1 = (a ^ 1 / 3)Dox0 = (x1)x1 = (x0 - (x0 ^ 3 - a) / (3 * x0 ^ 2))Loop While (Abs(x1 - x0)) > 0.000005 Print "a的立方根为："; x1End Sub3.编写牛顿法求方程根1)x3-x2-2x-3=0(初值x0=2)Private Sub Form_Click()Dim x0 as Single,x1 as Singlex1=2Dox0=x1x1=x0-(x0^3-x0^2-2*x0-3)/(3*x0^2-2*x0-2) Loop while abs(x1-x0)>0.00001Print x1End sub2)x-sinx=1/2(初值x0=1)Private Sub Form_Click()Dim x0 As Single, x1 As Singlex1 = 1Dox0 = x1x1 = x0 - (x0 - Sin(x0) - 1 / 2) / (1 - Cos(x0))Loop While Abs(x1 - x0) > 0.00001Print x1End Sub4.列主元高斯消去法Private Sub Form_Click()Dim a(1 To 3, 1 To 4) As Single, t#, i!, j!, k!, r!, l#, x(1 To 3) As Single For i = 1 To 3For j = 1 To 4a(i, j) = InputBox("输入一个数")Print a(i, j);Next jPrintNext iFor k = 1 To 2r = kFor i = k + l To 3If Abs(a(i, k)) > Abs(a(r, k)) Then r = iNext iIf r <> k ThenFor i = 1 To 4t = a(k, i)a(k, i) = a(r, i)a(r, i) = tNext iEnd IfFor i = k + l To 3l = (a(i, k) / a(k, k))For j = k + l To 4a(i, j) = (a(i, j) - l * a(k, j)) Next jNext iNext kFor k = 3 To 1 Step -1s = 0For j = k + l To 3s = s + (a(k, j) * x(j)) Next jx(k) = (a(k, 4) - s) / a(k, k) Next kFor i = 1 To 3Print x(i),Next iEnd sub5.LU分解法Private Sub Form_Click()Const n = 4Dim a(1 To n, 1 To n) As Single, l(1 To n, 1 To n) As Single, u(1 To n, 1 To n) As SingleDim x(1 To n) As Single, y(1 To n) As Single, b(1 To n) As Single, s#, i!, j!, k!, r!For i = 1 To nFor j = 1 To na(i, j) = InputBox("输入a数组")Print a(i, j)Next jPrintNext iFor i = 1 To nb(i) = InputBox("输入b数组")Print b(i)PrintFor k = 1 To nFor j = k To ns = 0For r = 1 To k - 1s = s + l(k, r) * u(r, j)Next ru(k, j) = a(k, j) - sNext jFor i = k + 1 To ns = 0For r = 1 To k - 1s = s + (l(i, r) * u(r, k)) Next rl(i, k) = (a(i, k) - s) / u(k, k) Next iNext kFor i = 1 To ns = 0For k = 1 To i - 1s = s + l(i, k) * y(k)y(i) = b(i) - sNext iFor i = n To 1 Step -1s = 0For k = i + 1 To ns = s + (u(i, k) * x(k))Next kx(i) = (y(i) - s) / u(i, i)Next iFor i = 1 To nPrint x(i)Next iEnd Sub6.雅克比迭代Option Base 1Function cha(x!(), y!()) As SingleDim z As Single, i As Single, k As Singlen = 3z = Abs(x(1) - y(1))For i = 2 To nIf (z < Abs(y(i) - x(i))) Then z = Abs(x(i) - y(i))Next icha = zEnd FunctionPrivate Sub Form_Click()Dim a1, x(3) As Single, y(3) As SingleDim t As Single, s As Single, a(3, 3) As SingleDim i As Integer, j As Integer, k As Integer, n As Integer n = 3a1 = Array(10, -2, -1, -2, 10, -1, -1, -2, 5)b = Array(3, 15, 10)For i = 1 To n: y(i) = 0: Next ik = 1For i = 1 To 3For j = 1 To 3a(i, j) = a1(k)k = k + 1Next j, iFor k = 1 To 30For i = 1 To nx(i) = y(i)Next iFor i = 1 To nt = 0For j = 1 To nIf (i <> j )Then t = t + a(i, j) * x(j) Next jy(i) = ((b(i) - t) / a(i, i))Next iIf (cha(x, y) < 0.00001 )Then Print k;For i = 1 To nPrint y(i)Next iExit ForEnd IfNext kIf k > 30 Then Print "发散"End Sub7. 高斯－赛德尔迭代Option Base 1Function cha(x!(), y!()) As SingleDim z As Single, i As Single, k As Singlen = 3z = Abs(x(1) - y(1))For i = 2 To nIf z < Abs(y(i) - x(i)) Then z = Abs(x(i) - y(i))Next icha = zEnd FunctionPrivate Sub Form_Click()Dim a1, x(3) As Single, y(3) As SingleDim t As Single, s As Single, a(3, 3) As SingleDim i As Integer, j As Integer, k As Integer, n As Integer n = 3a1 = Array(10, -2, -1, -2, 10, -1, -1, -2, 5)b = Array(3, 15, 10)For i = 1 To n: x(i) = 0: Next ik = 1For i = 1 To 3For j = 1 To 3a(i, j) = a1(k)k = k + 1Next j, iFor k = 1 To 30For i = 1 To ny(i) = x(i)Next iFor i = 1 To nt = 0For j = 1 To nIf i <> j Then t = t + a(i, j) * x(j) Next jx(i) = (b(i) - t) / a(i, i)Next iIf cha(x, y) < 0.00001 Then Print k;For i = 1 To nPrint x(i)Next iEnd IfNext kIf k > 30 Then Print "发散"End Sub8.拉格朗日插值多项式，求在ｔ＝３.５处的函数值的近似值，节点由ｘ，ｙ数组给出Private Sub Form_Click()Const n = 3Dim p#, s!Dim x, y As Variantx = Array(1, 2, 3, 4)y = Array(4, 5, 14, 37)t = InputBox("input t ")p = 0For k = 0 To ns = 1For i = 0 To nIf (i <> k) Thens = s * ((t - x(i)) / (x(k) - x(i)))Next ip = p + (y(k) * s)Next kPrint pEnd Sub9.牛顿基本插值公式Private Sub Form_Click()Const n = 4Dim x(n) As Single, y(n) As Single, t#, p#, s# For i = 0 To (n)x(i) = InputBox("input x" & Trim(Str(i)))y(i) = InputBox("input y" & Trim(Str(i))) Next it = InputBox("input t")For k = 1 To nFor i = n To k(-１)y(i) = (y(i) - y(i - 1)) / ((x(i) - x(i - k)))Next ip = y(0)h = 1For i = 1 To nh = (h * (t - x(i - 1)))p = p + (h * y(i))Next iPrint "p="; pEnd Sub10.拟合Private Sub Form_Click()Dim l#, m#, n#, i%, j%, k%, t1#Dim x As Variant, y As Variantn = 7m = 2x = Array(0, 1, 2, 3, 4, 5, 6, 7)y = Array(0, 5, 3, 2, 1, 2, 4, 7)ReDim a(0 To m, 0 To m + 1) As Single, t(n) As Single For i = 0 To mFor k = 1 To ns = s + x(k) ^ i * y(k) Next ka(i, m + 1) = sFor j = 0 To ms = 0For k = 1 To ns = s + x(k) ^ (i + j) Next ka(i, j) = sNext jNext iFor i = 0 To mFor j = o To m + 1 Print a(i, j),Next jPrintNext iFor k = 0 To mr = kFor i = k + 1 To mIf Abs(a(i, k)) > Abs(a(r, k)) Then r = iEnd IfNext iIf r <> k ThenFor i = 0 To m + 1t1 = a(k, i)a(k, i) = a(r, i)a(r, i) = t1Next iEnd Ifl = 1For i = k + 1 To ml = a(i, k) / a(k, k)For j = k + 1 To m + 1a(i, j) = a(i, j) - l * a(k, j)Next jNext iNext kFor k = m To step - 1s = 0For j = k + 1 To ms = s + a(k, j) * t(j)Next jt(k) = (a(k, m + 1) - s) / a(k, k) Next kPrint "y="; t(0)For i = 1 To mIf t(i) >= 0 Then Print "+" Print t(i); "*x^;i;"Next iEnd Sub11.SimpsonFunction f!(x!)f = x + x * Exp(x)End FunctionPrivate Sub Form_Click() Dim b!, N!, h!, x!, s!, a!a = 0:b = 1: N = 8h = (b - a) / (2 * N)x = as = f(a)For i = 1 To Nx = x + hs = s + 4 * f(x)x = x + hs = s + 2 * f(x)Next is = h / 3 * (s - f(b))Print sEnd Sub（2）Private Sub form_click()Dim a As Single, b As Single, eps As Single, s As Single Dim x As Single, h As Single, t1 As Single, t As Singlea = InputBox("输入积分下限a")b = InputBox("输入积分上限b")eps = InputBox("输入精度要求eps")h = b - at2 = (h / 2) * (f(a) + f(b))Dot1 = t2s = 0For x = a + h / 2 To b Step hs = s + f(x)Next xt2 = t1 / 2 + (h / 2) * sh = h / 2Loop While Abs(t1 - t2) > eps Print "积分的近似值："; t2End SubFunction f(x As Single) As Singlef = 1 / (1 + x * x)End Function12．用牛顿法求方程的附近根Private Sub Form_Click()Dim x#, x1#x1 = 1Dox = x1x1 = x - (x - Exp(-x)) / (1 + Exp(-x)) Loop While Abs(x1 - x) > 10 ^ (-5) Print x1End Sub。
( 鞋盒)上机纸张大小计算方法
（鞋盒）上机纸张大小计算方法一、3k版1）普通式：长=2底短折边+1底长+2底高+2盖长折边+1盖宽+2盖高宽=1底宽+2底高+2底长折边备注：（1底宽+2底高+2底长折边）≤（1盖长+2盖高+2盖短折边）时，则上述宽=1盖长+2盖高+2盖短折边：长=1底长+（1底宽+5cm）+ 1盖宽+2盖高+2盖长折边宽=2底长折边+1底宽+2底高备注：A、当2底高+2底短折边≥1底宽+5cm时，同1）算法。
B、当2底高+2底短折边≤1底宽+5cm时，同２）算法。
二、４k版1)普通式： 长＝１底长＋２底高＋２底短折边＋２盖短折边＋１盖长＋２盖高宽＝１底宽＋２底高＋２底长折边＋１盖宽＋２盖高＋２盖长折边２）加补强：长＝１底长＋２底高＋２底短折边＋２盖短折边＋１盖长＋２盖高 宽＝２底宽＋３底高＋３底长折边３）补强而且中间密合：长＝２底长＋２底高＋１底宽＋２底短折边 备注：（１底高＋１短折边）≤１底宽 宽＝２底宽＋３底高＋３底长折边４）交叉补强： 长＝２底长＋２底高＋（１底宽＋５cm）+２底短折边 宽＝２底宽＋３底高＋３底长折边备注：（当1底长折边+1底宽+1底高）≤（2盖高+1盖宽+2盖长折边）时，则2）、3）、4）三项宽=1底高+1底长折边+4盖高+4盖长折边+2盖宽。
三、６k版：长＝３底宽＋６底高＋６底长折边宽＝１底长＋２底短折边＋２底高＋１盖宽＋２盖高＋２盖长折边四、斜背式： １）３k 版：长＝１底宽＋１盖宽＋３底高＋１斜高＋４折边宽＝１盖长＋１底高＋１斜高＋２折边２）４k版：长＝２底长＋２底高＋２斜高＋４折边 宽＝１底宽＋１盖宽＋３底高＋１斜高＋４折边五、常规知识：１）底尺寸：４５０p、500p、550p、６００p；盖长＝底长+0.7cm;盖宽＝底宽＋０.5cm.坑纸则分别为＋0.8cm、＋０.6cm.２）盖尺寸：４５０p、500p、550p、６００p；底长＝盖长－0.7cm;底宽＝盖宽-０.5cm.坑纸则分别为-0.8cm、-０.6cm.３）斜背式内盒规定：斜高为：mens段为7.6cm,womens段为6.4cm,girls/boys段为5cm.折边为底高之四分之一。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
计算方法上机报告姓名：学号：班级：上课班级：说明：本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。
1. 对以下和式计算：∑∞⎪⎭⎫ ⎝⎛+-+-+-+=0681581482184161n n n n S n，要求：① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；（1） 算法思想1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 142111416818485861681n n n a n n n n n ε⎛⎫=---<< ⎪+++++⎝⎭； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)（2）算法结构1. ;0=s⎪⎭⎫⎝⎛+-+-+-+=681581482184161n n n n t n； 2. for 0,1,2,,n i =⋅⋅⋅ if 10m t -≤end;3. for ,1,2,,0n i i i =--⋅⋅⋅;s s t =+（3）Matlab源程序clear; %清除工作空间变量clc; %清除命令窗口命令m=input('请输入有效数字的位数m='); %输入有效数字的位数s=0;for n=0:50t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));if t<=10^(-m) %判断通项与精度的关系break;endend;fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值for i=n-1:-1:0t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6));s=s+t; %求和运算ends=vpa(s,m) %控制s的精度（4）结果与分析当保留11位有效数字时，需要将n值加到n=7，s =3.1415926536；当保留30位有效数字时，需要将n值加到n=22,s =3.14159265358979323846264338328。
通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。
2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。
在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。
已探测到一组等分点位置的深度数据(单位:米)如下表所示：①请用合适的曲线拟合所测数据点；②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；（1）算法思想如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。
分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。
在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。
海底光缆线的长度预测模型如下所示，光缆从A点铺h。
至B点，在某点处的深度为i海底光缆线的长度预测模型计算光缆长度时，用如下公式：20()L f x ds =⎰20'20()1()f x f x dx =+⎰ 191'20()1()k kk f x f x dx+==+∑⎰22()()x y =∆+∆（2）算法结构1. For n i ,,2,1,0⋅⋅⋅=1.1 i i M y ⇒2. For 2,1=k2.1 For k n n i ,,1, -=2.1.1 i k i i i i M x x M M ⇒----)/()(13. 101h x x ⇒-4. For 1-,,2,1n i =4.1 11++⇒-i i i h x x4.2 b a c c h h h i i i i i i ⇒⇒-⇒+++2;1;)/(114.3 i i d M ⇒+165. 0000;;c M d M d n n ⇒⇒⇒λn n n b a b ⇒⇒⇒2;;20μ6. 1111,γμ⇒⇒d b7. 获取M 的矩阵元素个数，存入m 8. For m k ,,3,2 =8.1 k k k l a ⇒-1/μ 8.2 k k k k c l b μ⇒⋅-1- 8.3 k k k k l d γγ⇒⋅-1- 9. m m m M ⇒μγ/10. For 1,,2,1 --=m m k10.1 k k k k k M M c ⇒⋅-+μγ/)(1 11. 获取x 的元素个数存入s 12. k ⇒113. For 1,,2,1-=s i13.1 if i x x ≤~then k i ⇒；break else k i ⇒+114. xx x x x x h x x k k k k ˆ~;~;11⇒-⇒-⇒--- y h x h M y x h M y x M x M k k k k k k ~/]ˆ)6()6(6ˆ6[2211331⇒-+-++---（3）Matlab 源程序clear; clc;x=0:1:20; %产生从0到20含21个等分点的数组X=0:0.2:20;y=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7. 95,8.86,9.81,10.80,10.93]; %等分点位置的深度数据n=length(x); %等分点的数目N=length(X);%% 求三次样条插值函数s(x)M=y;for k=2:3; %计算二阶差商并存放在M中for i=n:-1:k;M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));endendh(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;h(i)=x(i+1)-x(i);c(i)=h(i)/(h(i)+h(i-1));a(i)=1-c(i);b(i)=2;d(i)=6*M(i+1);endM(1)=0; %选择自然边界条件M(n)=0;b(1)=2;b(n)=2;c(1)=0;a(n)=0;d(1)=0;d(n)=0;u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);for k=2:n;l(k)=a(k)/u(k-1);u(k)=b(k)-l(k)*c(k-1);y1(k)=d(k)-l(k)*y1(k-1);endM(n)=y1(n)/u(n); %追赶法求解样条参数M(i) for k=n-1:-1:1;M(k)=(y1(k)-c(k)*M(k+1))/u(k);ends=zeros(1,N);for m=1:N;k=1;for i=2:n-1if X(m)<=x(i);k=i-1;break;elsek=i;endendH=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);x2=X(m)-x(k);s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/ 6))*x2)/H;end%% 计算所需光缆长度L=0; %计算所需光缆长度for i=2:NL=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);enddisp('所需光缆长度为L=');disp(L);figureplot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置','fontsize',16); %标注坐标轴含义ylabel('深度/m','fontsize',16);title('铺设河底光缆的曲线图','fontsize',16);grid;（4）结果与分析铺设海底光缆的曲线图如下图所示：仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为L=26.4844m。
3. 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。
时刻0 1 2 3 4 5 6 7 8 9 10 11 12平均气温15 14 14 14 14 15 16 18 20 20 23 25 28时刻13 14 15 16 17 18 19 20 21 22 23 24平均气温31 34 31 29 27 25 24 22 20 18 17 16（1）算法思想在本题中，数据点的数目较多。
当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。
用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。
另一方面，在有的实际问题中，用插值方法并不合适。
当数据点的数目很大时，要求)(x p 通过所有数据点，可能会失去原数据所表示的规律。
如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。
在这种情况下，不用插值标准而用其他近似标准更加合理。
通常情况下，是选取i α使2E 最小，这就是最小二乘近似问题。
在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2)(c t b ae C --=，计算相应的系数，估算误差，并作图比较各种函数之间的区别。
（2）算法结构本算法用正交化方法求数据的最小二乘近似。
假定数据以用来生成了G ，并将y 作为其最后一列（第1+n 列）存放。
结果在a 中，ρ是误差22E 。
I 、使用二次函数、三次函数、四次函数拟合时1. 将“时刻值”存入i x ，数据点的个数存入m2. 输入拟合多项式函数)(x p 的最高项次数1-=n h ，则拟合多项式函数为)()()()(2211x g x g x g x p n n ααα+⋅⋅⋅++= , 根据给定数据点确定G For1,,2,1,0-⋅⋅⋅=n jFor m i ,,2,1⋅⋅⋅=2.1 1,+⇒j i ji g x 2.2 1,+⇒n i i g y3. For n k ,,2,1⋅⋅⋅=3.1 [形成矩阵k Q ]3.1.1 σ⇒-∑=mk i ikkk g g 2/12))(sgn( 3.1.2 k kk g ωσ⇒-3.1.3 For m k k j ,,2,1⋅⋅⋅++=3.1.3.1 j jk g ω⇒ 3.1.4 βσω⇒k 3.2 [变换1-k G 到k G ] 3.2.1 kk g ⇒σFor 1,,,2,1+⋅⋅⋅++=n n k k j 3.2.2 t g mk i ij i ⇒∑=βω/)(3.2.3 For m k k i ,,1,⋅⋅⋅+=3.2.3.1 ij i ij g t g ⇒+ω4. [解三角方程1h Ra =] 4.1nnn n n a g g ⇒+/1.4.2 For 1,,2,1⋅⋅⋅--=n n i 4.2.1iii ni j j ijn i a g x gg ⇒-∑+=+/][11.5. [计算误差22E ]ρ⇒∑+=+mn i n i g121,II 、使用指数函数拟合时现将指数函数进行变形：将y C =,x t =代入2)(c t b aeC --=得:2)(c x b aey --=对上式左右取对数得: 222ln ln bx bcx bc a y -+-=令b bc bc a y z -==-==21202ln ln ααα，，，则可得多项式: 2210x x z ααα++=（3）Matlab 源程序clear; %清除工作空间变量 clc; %清除命令窗口命令 x=0:24; %将时刻值存入数组 y=[15,14,14,14,14,15,16,18,20,20,23,25,28,31,34,31,29,27,25,24,22,20,18,17,16]; [~,m]=size(x); %将数据点的个数存入m T=sum(y(1:m))/m;fprintf('一天的平均气温为T=%f\n',T); %求一天的平均气温 %% 二次、三次、四次函数的最小二乘近似h=input('请输入拟合多项式的最高项次数='); %根据给定数据点生成矩阵G n=h+1; G=[];for j=0:(n-1)g=x.^j; %g(x)按列排列 G=vertcat(G,g); %g 垂直连接G endG=G'; %转置得到矩阵Gfor i=1:m %将数据y 作为G 的最后一列（n+1列）G(i,n+1)=y(i);endG;for k=1:n %形成矩阵Q(k) if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gk for j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA (n)=G(n,n+1)/G(n,n); %解三角方程求系数Afor i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);ende=sum(G(n+1:m,n+1).^2); %计算误差efprintf('%d次函数的系数是：',h); %输出系数a及误差e disp(A);fprintf('使用%d次函数拟合的误差是%f：',h,e);t=0:0.05:24;A=fliplr(A); %将系数数组左右翻转Y=poly2sym(A); %将系数数组转化为多项式subs(Y,'x',t);Y=double(ans);figure(1)plot(x,y,'k*',t,Y,'r-'); %绘制拟合多项式函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title([num2str(n-1),'次函数的最小二乘曲线']);grid;%% 指数函数的最小二乘近似yy=log(y);n=3;G=[];GG=[];for j=0:(n-1)g=x.^j; %g(x)按列排列G=vertcat(G,g); %g垂直连接Ggg=t.^j; %g(x)按列排列GG=vertcat(GG,gg); %g垂直连接GendG=G'; %转置得到矩阵Gfor i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=yy(i);endG;for k=1:n %形成矩阵Q(k)if G(k,k)>0;sgn=1;elseif G(k,k)==0;sgn=0;else sgn=-1;endsgm=-sgn*sqrt(sum(G(k:m,k).^2));w=zeros(1,n);w(k)=G(k,k)-sgm;for j=k+1:mw(j)=G(j,k);endbt=sgm*w(k);G(k,k)=sgm; %变换Gk-1到Gkfor j=k+1:n+1t=sum(w(k:m)*G(k:m,j))/bt;for i=k:m;G(i,j)=G(i,j)+t*w(i);endendendA(n)=G(n,n+1)/G(n,n); %解三角方程求系数A for i=n-1:-1:1A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);endb=-A(3);c=A(2)/(2*b);a=exp(A(1)+b*(c^2));G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));e=sum(G(n+1:m,n+1).^2); %计算误差efprintf('\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c); %输出系数及误差e fprintf('\n使用指数函数拟合的误差是：%f',e);t=0:0.05:24;YY=a.*exp(-b.*(t-c).^2);figure(2)plot(x,y,'k*',t,YY,'r-'); %绘制拟合指数函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');title(['指数函数的最小二乘曲线']);grid;（4）结果与分析a、二次函数：一天的平均气温为：21.20002次函数的系数: 8.3063 2.6064 -0.0938使用2次函数拟合的误差是：280.339547二次函数的最小二乘曲线如下图所示：b、三次函数：一天的平均气温为：21.20003次函数的系数: 13.3880 -0.2273 0.2075 -0.0084 使用3次函数拟合的误差是：131.061822三次函数的最小二乘曲线如下图所示：c、四次函数：一天的平均气温为：21.20004次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009 使用4次函数拟合的误差是：59.04118四次函数的最小二乘曲线如下图所示：d、指数函数：一天的平均气温为：21.2000指数函数的系数是：a=26.160286，b=0.004442，c=14.081900使用指数函数拟合的误差是：57.034644指数函数的最小二乘曲线如下图所示：通过上述几种拟合可以发现，多项式的次数越高，计算拟合的效果越好，误差越小，说明结果越准确；同时，指数多项式拟合的次数虽然不高，但误差最小，说明结果最准确。