计算方法上机作业

合集下载

西安交通大学计算方法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

计算方法与计算实验一误差分析

（1）MATLAB 主程序 function [k,juecha,xiangcha,xk]= liti112(x0,x1,limax) % 输入的量--x0是初值, limax是迭代次数和精确值x;
% 输出的量--每次迭代次数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.折边为底高之四分之一。

《计算方法》上机实验指导书刘轶中-8页word资料

理学院《计算方法》实验指导书适合专业：信息与计算科学数学与应用数学贵州大学二OO七年八月前言《计算机数值计算方法》包括很多常用的近似计算的处理手段和算法，是计算科学与技术专业的必修课程，为了加强学生对该门课程的理解，使学生更好地掌握书中的计算方法、编制程序的能力，学习计算方法课程必须重视实验环节，即独立编写出程序，独立上机调试程序，必须保证有足够的上机实验时间。

在多年教学实践基础上编写了《计算机数值计算方法》上机实习指南，目的是通过上机实践，使学生能对教学内容加深理解，同时培养学生动手的能力.本实习指南，可与《计算机数值计算方法》课本配套使用，但是又有独立性，它不具体依赖哪本教科书，主要的计算方法在本指南中都有，因此，凡学习计算方法课的学生都可以参考本指南进行上机实习。

上机结束后，按要求整理出实验报告。

实验报告的内容参阅《计算机数值计算方法》上机实验大纲。

目录第一章解线性方程组的直接法实验一 Gauss列主元素消去法实验二解三对角线性方程组的追赶法第二章插值法与最小二乘法实验三 lagrange插值法实验四分段插值法实验五曲线拟合的最小二乘法第三章数值积分实验六复合求积法实验七变步长法第四章常微分方程数值解法实验八 Euler 方法第五章解线性方程组和非线性方程的迭代法实验九 Jacobi 迭代法、Gauss-Seidel 迭代法实验十 Newton 迭代法实验一： Gauss 列主元素消去法实验学时：2实验类型：验证实验要求：必修一、实验目的用gauss 消去法求线性方程组AX=b. 其中一、实验内容二、实验条件PC 机，tc2.0，Internet 网。

三、实验步骤1．根据算法事先写出相应程序。

2．启动PC 机，进入tc 集成环境，输入代码。

3．编译调试。

4．调试通过，计算出正确结果后。

实验二解三对角线性方程组的追赶法⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=b b b x x x a a a a a a a a a n n nn n n n n b X A M M 2122122221112111....................................实验学时：2实验类型：验证实验要求：必修一、实验目的二、实验内容三、实验组织远行要求统一进行实验，一人一组四、实验条件PC机，tc2.0，Internet网五、实验步骤a)根据算法事先写出相应程序。

计算方法上机作业

方程求根一、目的和意义非线性方程在科学研究与工程实践中广泛出现，例如，优化问题、特征值问题、微分方程问题等。

但是，除少量方程外，大多数非线性方程求根相当困难，常见的几个简单、有效的数值求根方法，包括二分法、迭代法、牛顿法、割线法等，本实验旨在比较各种算法的计算性能和使用范围。

二、计算公式1.二分法2.不动点迭代法三、结构程序设计代码1.二分法1).定义所求解函数function [ y ] = f( x )y = x^3 + 4*x^2 - 10;end2).执行算法%初始化，设置区间端点a、b，误差限tola = 1;b = 2; tol = 0.5*10^(-6);k = 0; fa = f(a);%设置最大二分次数为30for k = 1:50p = (a + b)/2; fp = f(p);if(fp == 0 || (b - a)/2 < tol)breakendif(fa * fp < 0)b = p;elsea = p;enddisp('近似解p = ');disp(vpa(p,10)); disp('迭代次数k = ');disp(k);end2.不动点迭代法1).定义不动点方程g(x)function [ y ] = g( x )y = x^3 + 4*x^2 + x - 10;end2).执行算法%初始化，设置误差限，设置初值p0tol = 0.5*10^(-6);k = 1; p0 = 1.5;%迭代次数为10次while k <= 10p = g(p0);if abs(p - p0) < tolbreakenddisp('近似解p = ');disp(vpa(p,10)); disp('迭代次数k = ');disp(k);k = k + 1; p0 = p;end四、结果及其讨论1.二分法结果由于结果较长，只取了一部分，从图中可以看出，迭代20次可得到误差限范围内的近似解p=1.36522。

数值计算方法上机实验报告

数值计算方法上机实验报告实验目的：复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。

利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。

上机练习任务：利用计算机基本C 语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。

掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。

一、各算法的算法原理及计算机程序框图1. 列主元高斯消去法算法原理：高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘一个方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上对上三角方程组求解。

列选住院是当高斯消元到第k 步时，从k 列的kk a 以下（包括kk a ）的各元素中选出绝对值最大的，然后通过行交换将其交换到kk a 的位置上。

交换系数矩阵中的两行（包括常数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结果。

●源程序：#define N 200#include "stdio.h"#include "math.h"FILE *fp1,*fp2;void LZ(){int n,i,j,k=0,l;double d,t,t1;static double x[N],a[N][N];fp1=fopen("a1.txt","r");fp2=fopen("b1.txt","w");fscanf(fp1,"%d",&n);for(i=0;i<n;++i)for(j=0;j<=n;++j){fscanf(fp1,"%lf",&a[i][j]);}{d=a[k][k];l=k;i=k+1;do{if(fabs(a[i][k])>fabs(d)) /*选主元*/{d=a[i][k];l=i;}i++;}while(i<n);if(d==0){printf("\n输入矩阵有误！\n");}else{ /*换行*/if(l!=k){for(j=k;j<=n;j++){t=a[l][j];a[l][j]=a[k][j];a[k][j]=t;}}}for(j=k+1;j<=n;j++) /*正消*/ a[k][j]/=a[k][k];for(i=k+1;i<n;i++)for(j=k+1;j<=n;j++)a[i][j]-=a[i][k]*a[k][j];k++;}while(k<n);if(k!=0){for(i=n-1;i>=0;i--) /*回代*/ {t1=0;for(j=i+1;j<n;j++)t1+=a[i][j]*x[j];x[i]=a[i][n]-t1;}for(i=0;i<n;i++)fprintf(fp2,"\n 方程组的根为x[%d]=%lf",i+1,x[i]); fclose(fp1); fclose(fp2); }main() { LZ(); }● 具体算例及求解结果：用列选主元法求解下列线性方程组⎪⎩⎪⎨⎧=++=++=-+28x x 23x 2232832321321321x x x x x x 输入3 输出结果：方程组的根为x[1]=6.0000001 2 -3 8 方程组的根为x[2]=4.000000 2 1 3 22 方程组的根为x[3]=2.000000 3 2 1 28● 输入变量、输出变量说明：输入变量：ij a 系数矩阵元素，i b 常向量元素输出变量：12,,n b b b 解向量元素2. 杜里特尔分解法解线性方程● 算法原理：求解线性方程组Ax b =时，当对A 进行杜里特尔分解，则等价于求解LUx b =，这时可归结为利用递推计算相继求解两个三角形（系数矩阵为三角矩阵）方程组，用顺代，由Ly b =求出y ，再利用回带，由Ux y =求出x 。

数值计算方法上机题

习题二问题：1．编制通用子程序对n+1个节点x i及y i=f(x i) (i=1,…n) （1）n次拉格朗日插值计算公式；（2）n次牛顿向前插值计算公式；（3）n次牛顿向后插值计算公式；（一）程序流程图（1）拉格朗日插值程序流程图（2）牛顿向前插值程序流程图（3）牛顿向后插值程序流程图。

（二）源程序见主程序清单问题：2．计算（1）已知f(x)=lnx,，[a，b]=[1，2]，取h=0.1，x i=1+ih，i=0，1， (10)用通用程序（1），（3）计算ln1.54及ln1.98的近似值；（一）程序清单/* program of question 2.1, page 61 */#include "stdio.h"#include "math.h"main(){ int i,flag=0;double z1,z2,x[11],y[11],t,s1,s2,z[10],c[11][11],log(double),ntb(),L(); for(i=0;i<=10;i++){x[i]=1+0.1*i;y[i]=log(x[i]);}printf("data x:\n");for(i=0;i<=10;i++){flag++;printf("%11.6f",x[i]);if(flag%4==0)printf("\n");}printf("\ndata y:\n");flag=0;for(i=0;i<=10;i++){flag++;printf("%11.6f",y[i]);if(flag%4==0)printf("\n");}printf("\nThe true value:\n");printf(" ln1.54=%f ln1.98=%f\n",log(1.54),log(1.98));z1=L(x,y,10,1.54);z2=L(x,y,10,1.98);t=(1.54-x[10])/0.1;s1=ntb(y,10,t,z,c);s2=ntb(y,10,t,z,c);t=(1.98-x[10])/0.1;printf("The approximate value:\n");printf(" L(1.54)=%f L(1.98)=%f\n",z1,z2);printf(" NTB(1.54)=%f NTB(1.98)=%f\n",s1,s2);}double L(double x[],double y[],int n,double t){ int i,k; double z=0.0,s;if(n==1)z=y[0];for(k=0;k<=n;k++){ s=1.0;for(i=0;i<=n;i++) if(i!=k)s=s*(t-x[i])/(x[k]-x[i]);z=z+s*y[k]; }return z;}double ntb(double y[],int n,double t,double z[],double c[][11]) { int i,j,sn=n;double s;z[0]=t;for(i=1;i<=n-1;i++) z[i]=z[i-1]*(t+i)/(i+1);for(i=0;i<=n;i++) c[i][0]=y[sn--];for(j=1;j<=n;j++)for(i=0;i<=n-j;i++) c[i][j]=c[i][j-1]-c[i+1][j-1];s=y[n];for(i=0;i<=n-1;i++) s=s+z[i]*c[0][i+1];return s;}(二)运行结果data x:1.000000 1.100000 1.200000 1.3000001.400000 1.500000 1.600000 1.7000001.800000 1.9000002.000000data y:0.000000 0.095310 0.182322 0.2623640.336472 0.405465 0.470004 0.5306280.587787 0.641854 0.693147The true value:ln1.54=0.431782 ln1.98=0.683097The approximate value:L(1.54)=0.431782 L(1.98)=0.683097NTB(1.54)=0.431782 NTB(1.98)=0.683097问题：（2）f(x)=1/(1+25x2), |x|≤1取等距节点n=5和n=10，用通用程序（1），（2）依次计算x=-0.95+ih（i=0，1， (19)h=0.1）处f(x)的近似值，并将其结果与其真实值相比较。

数值计算方法上机实验题

数值计算方法上机实验实验内容：1.要求:分别用复化梯形法,复化Simpson 法和 Romberg 公式计算.2．给定积分dx e x⎰31和dx x ⎰311 ,分别用下列方法计算积分值要求准确到510- ,并比较分析计算时间. 1）变步长梯形法; 2）变步长 Simpson 法; 3) Romberg 方法.算法描述：1、复合梯形法：⎰=tdt t a t V 0)()( ))()(2)((211∑-=++=n k k n b f x f a f hT输入被积函数数据点t,a. 输出积分值.n T复合Simpson 法：⎰=tdt t a t V 0)()( ))()(2)(4)((6101121∑∑---=++++=n k n k k k n b f x f x f a f hS输入被积函数f(x),积分区间[a,b]和n 输出复合Simpson 积分值n S步1 .);()(;a x b f a f S nab h n ⇐-⇐-⇐ 步2 对n k ,,2,1 =执行).(2;2);(4;2x f S S hx x x f S S h x x n n n n +⇐+⇐+⇐+⇐步3 n n S hS ⨯⇐6步4 输出n SRomberg 积分法：根据已知数据对其进行多项式拟合得出p(x);f(x)⇐p(x); 输入被积函数f(x)，积分区间端点a,b,允许误差ε 输出 Romberg 积分值n R 2 步1 .0;0;0;0));()((2;1111⇐===+⇐-⇐k R C S b f a f hT a b h 步2 反复执行步3→步9. 步3 .2;0h a x S +⇐⇐ 步4 反复执行步5→步6. 步5 ;);(h x x x f S S +⇐+⇐步6 若x ≥b,则退出本层循环. 步7 执行.6316364;1511516;3134;2212212212212C C R S S C T T S S h T T -⇐-⇐-⇐+⇐步8 执行.1;;;;;2;2121212112+⇐⇐⇐⇐⇐⇐-⇐k k R R C C S S T T hh R R e 步9 若e ≤ε且k ≥5,则退出循环. 步10 .22R R n ⇐ 步11 输出.2n R2、变步长梯形算法：功能求积分⎰ba)(dx x f ，允许误差为ε。

计算方法大作业——三次样条插值

8
计算方法上机报告
此完成所有数据的输入。继续按 Enter 键会出现提示“选择封闭方程组的边界条件: 第一类边界条件输入 1，第二类边界条件输入 2，第三类边界条件输入 3。 ”根据已知情况选择相应的边界条件，若为自然三次样条插值，则选 1，并将插值区间两端点的二阶导数值设置为 0。输入完成之后按 Enter 开始求解，程序运行结束后命令窗口会显示要求的三次样条插值函数，同时会出现该插值函数以及插值节点的图像，便于直接观察。 2.3 算例及计算结果 (1) 《数值分析》课本第 137 页的例题 4.6.1，已知函数 y=f(x)的数值如下表，求它的自然三次样条插值函数。 xi yi -3 7 -1 11 0 26 3 56 4 29
(2) 给定函数 f ( x)
3 x 1 1 x 0 0 x3 3 x 4
1 (1 x 1) 。取等距节点，构造牛顿插值多项式 N5(x) 1 25x 2 和 N10(x)及三次样条插值函数 S10(x)。分别将三种插值多项式与 f(x)的曲线画在同一个
N10 x
22757 10 5444 8 20216 6 17147 4 3725 2 x x x x x 1 103 11 53 139 221
将牛顿插值多项式 N5(x)和 N10(x)及三次样条插值函数 S10(x)分别与 f(x)的曲线画在同一个坐标系上进行比较，如图 12。可以看出三次样条函数与原函数符合的非常好，对于低次的牛顿插值多项式，与原函数的大致趋势相同，而高次的牛顿插值多项式由于龙格现象的出现，与原函数之间相差比较大。
S ( xi ) S ( xi ), ( xi ) S ( xi ), S S ( x ) S ( x ), i i i 1, 2, , n 1

数值计算方法上机实习题答案.doc

1．设I n 1 x ndx ，0 5 x（ 1）由递推公式 I n 5I n 11，从 I 0的几个近似值出发，计算I 20；n解：易得： I 0 ln6-ln5=0.1823, 程序为：I=0.182;for n=1:20I=(-5)*I+1/n;endI输出结果为： I 20= -3.0666e+010（ 2）粗糙估计 I 20，用 I n 1 1I n 1 1 ，计算 I 0；5 5n0.0079 1 x 20 1 x 200.0095因为dx I 20dx 6 5所以取 I 20 1(0.0079 0.0095) 0.0087 2程序为： I=0.0087;for n=1:20I=(-1/5)*I+1/(5*n);endII 0= 0.0083（ 3）分析结果的可靠性及产生此现象的原因(重点分析原因 )。

首先分析两种递推式的误差；设第一递推式中开始时的误差为E0 I 0 I 0，递推过程的舍入误差不计。

并记 E n I n I n，则有 E n 5E n 1 ( 5) n E0。

因为 E20( 5) 20 E0 I 20，所此递推式不可靠。

而在第二种递推式中E0 1E1 (1)n E n，误差在缩小，5 5所以此递推式是可靠的。

出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制，即算法是否数值稳定。

2．求方程e x10x 2 0 的近似根，要求x k 1x k 5 10 4，并比较计算量。

（1）在 [0， 1]上用二分法；程序： a=0;b=1.0;while abs(b-a)>5*1e-4c=(b+a)/2;if exp(c)+10*c-2>0b=c;else a=c;endendc结果： c =0.0903（ 2）取初值x0 0，并用迭代 x k 1 2 e x ；10程序： x=0;a=1;while abs(x-a)>5*1e-4a=x;x=(2-exp(x))/10;endx结果： x =0.0905（3）加速迭代的结果；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;y=exp(x)+10*x-2;z=exp(y)+10*y-2;x=x-(y-x)^2/(z-2*y+x);b=x;endx结果： x =0.0995（ 4）取初值x00 ，并用牛顿迭代法；程序： x=0;a=0;b=1;while abs(b-a)>5*1e-4a=x;x=x-(exp(x)+10*x-2)/(exp(x)+10); b=x;end x 结果： x =0.0905（ 5）分析绝对误差。

中国石油大学(华东)计算方法上机实验程序

1.二分法#include <iostream.h>double f(double x);void main(){double a,b,x,ep;a=0;b=5;ep=0.000001;while(b-a>=ep){x=(a+b)/2;if(f(a)*f(x)>0) a=x;else b=x;}x=(a+b)/2;cout<<"the root of f(x) is"<<" "<<x<<endl; }double f(double x){double f1;f1=x*x*x-6*x-1;return f1;}2.牛顿迭代法#include <iostream.h>#include <math.h>double f1(double x);double f2(double x);void main(){int flag,nmax,k;double x,ep,x0;x=2;ep=000001;nmax=200;flag=1;k=1;while(fabs(x-x0)>ep&&(k<nmax)){x0=x;if(fabs(f2(x0))<ep){flag=0;break;}x=x0-f1(x0)/f2(x0);k=k+1;}if(k>=nmax) flag=0;if(flag==0)cout<<"the newton method is failure"<<endl; elsecout<<"the root of f(x) is"<<" "<<x<<endl;}double f1(double x){double f1;f1=x*x*x-3*x-1;return f1;}double f2(double x){double f1;f1=3*x*x-3;return f1;}3.列主元消去法#include<iostream.h>#include<math.h>void main(){static double a[3][3]={{7,8,11},{5,1,-3},{1,2,3}};static double b[3]={-3,-4,1};double ep=0.001;int ip=1;int n=3;double dmax,temp,s;int m,i,j,k;for(k=0;k<=n-2;k++){dmax=fabs(a[k][k]);m=k;for(i=k+1;i<=n-1;i++){if(fabs(a[i][k])>dmax){dmax=fabs(a[i][k]);m=i;}}if(dmax<ep){ip=-1;break;}if(m!=k){for(j=k;j<=n-1;j++){temp=a[k][j];a[k][j]=a[m][j];a[m][j]=temp;}temp=b[k];b[k]=b[m];b[m]=temp;}for(i=k+1;i<=n-1;i++){a[i][k]=a[i][k]/a[k][k];for(j=k+1;j<=n-1;j++){a[i][j]=a[i][j]-a[i][k]*a[k][j];}b[i]=b[i]-a[i][k]*b[k];}}b[n-1]=b[n-1]/a[n-1][n-1];for(i=n-2;i>=0;i--){s=0;for(j=i+1;j<=n-1;j++){s=s+a[i][j]*b[j];}b[i]=(b[i]-s)/a[i][i];}if(ip==-1)cout<<"gauss method is failure"<<endl;elsecout<<"the solution of eqution is "<<b[0]<<","<<b[1]<<","<<b[2]<<endl; }4.LU分解法#include <iostream>#include <fstream>using namespace std;#define N 3int main(){float A[N][N],U[N][N],L[N][N],X[N],Y[N],B[N];float s=0;int i,j,k;cout<<"请输入将要分解的矩阵:"<<endl;for(i=0;i<N;i++) //输入待分解的矩阵for(j=0;j<N;j++)cin>>A[i][j];for(i=0;i<N;i++) //初始化for(j=0;j<N;j++){L[i][j]=0;U[i][j]=0;if(i==j)L[i][j]=1;}for(j=0;j<N;j++) //求U矩阵第一行U[0][j]=A[0][j];for(i=1;i<N;i++)L[i][0]=A[i][0]/U[0][0]; //求L矩阵的第一列for(k=1;k<N;k++) //用k记录化解的行列数，在行和列分解都完成时再进行下一行列的分解{for(j=k;j<N;j++) //本循环将完成行分解的工作{for(int u=0;u<=k-1;u++){s=s+L[k][u]*U[u][j];}U[k][j]=A[k][j]-s;s=0; //每行的和不一致，一行使用完后要置0}for(i=k+1;i<N;i++) //在行分解完成后进行列分解{for(int l=0;l<=k-1;l++){s=s+L[i][l]*U[l][k];}L[i][k]=(A[i][k]-s)/U[k][k];s=0; //每列的和不一致，一行使用完后要置0 }}cout<<endl<<"L"<<endl;for(i=0;i<N;i++){ //输出L Ufor(j=0;j<N;j++)cout<<L[i][j]<<"\t";cout<<endl;}cout<<endl<<"U"<<endl;for(i=0;i<N;i++){for(j=0;j<N;j++)cout<<U[i][j]<<"\t";cout<<endl;}cout<<"请输入矩阵B"<<endl;cout<<endl;for(i=0;i<N;i++)cin>>B[i];Y[0]=B[0]; //求Yfor(k=1;k<N;k++){for(j=0;j<k;j++){s=s+L[k][j]*Y[j];}Y[k]=B[k]-s;s=0;}cout<<endl<<"Y:"<<endl; //输出Yfor(i=0;i<N;i++)cout<<Y[i]<<"\t";cout<<endl;X[N-1]=Y[N-1]/U[N-1][N-1]; //求Xfor(k=N-2;k>=0;k--){for(j=k+1;j<N;j++){s=s+U[k][j]*X[j];}X[k]=(Y[k]-s)/U[k][k];s=0;}cout<<endl<<"X:"<<endl; //输出Xfor(i=0;i<N;i++)cout<<X[i]<<"\t";cout<<endl;system("pause");return 0;}5.追赶法#include <iostream>using namespace std;int main(){double a[15],b[15],c[15],d[15];double t;int i,n;cout<<"请输入 n= ";cin>>n;cout<<endl;cout<<"请输入 b[1],c[1],d[1]:"<<endl;cin>>b[1]>>c[1]>>d[1];cout<<endl;cout<<"输入2行到n-1行 :"<<endl;for(i=2;i<=n-1;i++){cin>>a[i]>>b[i]>>c[i]>>d[i];}cout<<endl;cout<<"请输入 a[n],b[n],d[n]:"<<endl;cin>>a[n]>>b[n]>>d[n];c[1]=c[1]/b[1];d[1]=d[1]/b[1];for(i=2;i<=n-1;i++){t=b[i]-c[i-1]*a[i];c[i]=c[i]/t;d[i]=(d[i]-d[i-1]*a[i])/t;}d[n]=(d[n]-d[n-1]*a[n])/(b[n]-c[n-1]*a[n]);for(i=n-1;i>=1;i--) d[i]=d[i]-c[i]*d[i+1];for(i=1;i<=n;i++)cout<<endl<<"x["<<i<<"]="<<d[i];cout<<endl;return 0;}6.平方根法#include <stdio.h>#include <stdlib.h>#include <math.h>#define EPS 1.0e-8#define N 20double a[N][N], b[N], x[N];int n;int zhuyuan(int row); /* 选主元*/void hangjiaohuan(int row1, int row2); /* 行交换*/ void xiaoyuan(int row); /*消元*/void huidai(); /* 回代*/void main(){printf("请输入方程的维数n!\n n = ");scanf("%d", &n);getchar();printf("\n输入%d行%d列矩阵\n", n, n);for (int i=0; i<n; i++){for (int j=0; j<n; j++)scanf("%lf", &a[i][j]);getchar();}printf("\n输入线性方程组右端项b[%d]: ", n); for (i=0; i<n; i++){scanf("%lf", &b[i]);}getchar();for (i=0; i<n-1; i++){double rowmark = zhuyuan(i); if (rowmark == -1){printf("多解!");system("pause");return;}if (rowmark != i){hangjiaohuan(i, rowmark);}xiaoyuan(i);}huidai();printf("\n线性方程组的解为: "); for (i=0; i<n; i++){printf("x%d=%lf ", i+1, x[i]); }printf("\n");system("pause");}int zhuyuan(int row){double elem = a[row][row];int rowmark = row;for (int i=row+1; i<n; i++){if (elem<a[i][row]){elem = a[i][row];rowmark = i;}}if(fabs(elem) <= EPS){return -1;return rowmark;}void hangjiaohuan(int row1, int row2) {double tmp;tmp = b[row1];b[row1] = b[row2];b[row2] = tmp;for (int j=0; j<n; j++){tmp = a[row1][j];a[row1][j] = a[row2][j];a[row2][j] = tmp;}}void xiaoyuan(int row){for (int i=row+1; i<n; i++){int j=row;double tmp = a[i][j]/a[row][row]; b[i] -= tmp*b[row];for (a[i][j++] = 0; j<n; j++){a[i][j] -= tmp*a[row][j];}}}void huidai(){x[n-1] = b[n-1]/a[n-1][n-1];for (int i=n-2; i>=0; i--){double sum = 0.0;for (int j=i+1; j<n; j++){sum -= a[i][j]*x[j];}x[i] = (b[i]+sum)/a[i][i];}7.G-S迭代法#include<iostream>#include<math.h>void main(){double ep=0.000001;int n=3,ip=1,nmax=200;double c[3][3]={{0,0.2,0.1},{0.2,0,0.1},{0.2,0.4,0}};static double d[3]={0.3,1.5,2};static double x[3]={0.3,1.5,2};int i,j,k;double emax,s;k=0;do{emax=0;for(i=0;i<n;i++){s=d[i];for(j=0;j<n;j++){s=s+c[i][j]*x[j];}if(fabs(s-x[i])>emax)emax=fabs(s-x[i]);x[i]=s;}k=k+1;}while((emax>ep)&&(k<nmax));if(k>=nmax)ip=-1;elseip=1;if(ip==1)printf("the soulation of eqution is %1.10f,%1.10f,%1.10f",x[0],x[1],x[2]);else printf("failure");}8.曲线拟合法#include <iostream>#include <cmath>#define gasvoid main(){static double x[10]={1.36,1.49,1.73,1.81,1.95,2.16,2.28,2.48};static double y[10]={14.094,15.069,16.844,17.378,18.435,19.949,20.963,22.495};int n=8;int m=2;double a[10][10],b[10],z[30];int i,j,k,ip;double s;z[0]=n;for(i=0;i<=2*m-1;i++){s=0;for(k=0;k<=n-1;k++){s=s+pow(x[k],(i+1));}z[i+1]=s;}for(i=0;i<=m;i++){s=0;for(k=0;k<=n-1;k++){s=s+y[k]*pow(x[k],(i));}b[i]=s;}for(i=0;i<=m;i++){for(j=0;j<=m;j++){a[i][j]=z[i+j];}}gas(a,b,m+1,ip,0.000001);if(ip!=-1)printf("%f+%fx+%fx*x",b[0],b[1],b[2]);elseprintf("failure"); }。

计算方法实验上机报告(完整版)

简单迭代法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(12*x+sin(x)-1,1.0/3)void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果牛顿迭代法一#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)-sin(x)-12*x+1)/(3*pow(x,2)-cos(x)-12) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果埃特金加速法#include<stdio.h>#include<math.h>#define x0 3.0#define MAXREPT 1000#define EPS 1E-6#define G(x) (pow(x,3)-sin(x)+1)/12void main(){int i;double x1=x0,x2=x0;double z,y;printf("k\tx1\tx2\txk\n");for(i=0;i<MAXREPT;i++){if(i==0)printf("%d\t\t\t%g\n",i,x2);elseprintf("%d\t%g\t%g\t%g\n",i,y,z,x2);y=G(x1);z=G(y);x2=z-((z-y)*(z-y))/(z-2*y+x1);if (fabs(x2-x1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x2,i);return;}x1=x2;}printf("AFTER %d repeate,no solved.\n",MAXREPT);} 结果牛顿迭代法二#include<stdio.h>#include<math.h>#define x0 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) x-(pow(x,3)+pow(x,2)-3*x-3)/(3*pow(x,2)+2*x-3) void main(){int i;double x_k=x0,x_k1=x0;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k1);x_k1=G(x_k);if (fabs(x_k1-x_k)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k1,i);return;}x_k=x_k1;}printf("AFTER %d repeate,no solved.\n",MAXREPT);}结果弦截法#include<stdio.h>#include<math.h>#define x0 0#define x1 1.5#define MAXREPT 1000#define EPS 1E-6#define G(x) pow(x,3)+pow(x,2)-3*x-3void main(){int i;double x_k=x0,x_k1=x1,x_k2=0;double y,z;printf("k\txk\n");for(i=0;i<MAXREPT;i++){printf("%d\t%g\n",i,x_k2);y=G(x_k);z=G(x_k1);x_k2=x_k1-(z*(x_k1-x_k))/(z-y);if (fabs(x_k2-x_k1)<EPS){printf("THE ROOT IS x=%g,k=%d\n",x_k2,i);return;}x_k=x_k1;x_k1=x_k2;}printf("AFTER %d repeate,no solved.\n",MAXREPT); } 结果高斯顺序消元法#include<stdio.h>#include<math.h>#define N 4static double aa[N][N+1]={{2,4,0,1,1},{3,8,2,2,3},{1,3,3,0,6},{2,5,2,2,3}}; int gauss(double a[][N+2],double x[]);void putout(double a[][N+2]);void main(){int i,j,det;double a[N+1][N+2],x[N+1];for(i=1;i<=N;i++)for(j=1;j<=N+1;j++)a[i][j]=aa[i-1][j-1];det=gauss(a,x);if(det!=0)for(i=1;i<=N;i++)printf(" x[%d]=%g",i,x[i]);printf("\n");}int gauss(double a[][N+2],double x[]){int i,j,k;double c;putout(a);for(k=1;k<=N-1;k++){ if(fabs(a[k][k])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(i=k+1;i<=N;i++){c=a[i][k]/a[k][k];for(j=k;j<=N+1;j++){a[i][j]=a[i][j]-c*a[k][j];}}putout(a);}if(fabs(a[N][N])<1e-17){printf("\n pivot element is 0.fail!\n");return 0;}for(k=N;k>=1;k--){x[k]=a[k][N+1];for(j=k+1;j<=N;j++){x[k]=x[k]-a[k][j]*x[j];}x[k]=x[k]/a[k][k];}return(1);}void putout(double a[][N+2]){for(int i=1;i<=N;i++){for(int j=1;j<=N+1;j++)printf("%-15g",a[i][j]);printf("\n");}printf("\n");}结果雅克比迭代法#include<stdio.h>#include<math.h>#define N 5#define EPS 0.5e-4static double aa[N][N]={{4,-1,0,-1,0},{-1,4,-1,0,-1},{0,-1,4,-1,0},{-1,0,-1,4,-1},{0,-1,0,-1,4}}; static double bb[N]={2,1,2,1,2};void main(){int i,j,k,NO;double a[N+1][N+1],b[N+1],x[N+1],y[N+1];double d,sum,max;for(i=1;i<=N;i++){for(j=1;j<=N;j++)a[i][j]=aa[i-1][j-1];b[i]=bb[i-1];}printf("\n 请输入最大迭代次数（尽量取大值！）NO:");scanf("%d",&NO);printf("\n");for(i=1;i<=N;i++)x[i]=0;k=0;printf(" k",' ');for(i=1;i<=N;i++)printf("%8cx[%d]",' ',i);printf("\n 0");for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");do{for(i=1;i<=N;i++){sum=0.0;for(j=1;j<=N;j++)if(j!=i) sum=sum+a[i][j]*x[j];y[i]=(-sum+b[i])/a[i][i];}max=0.0;for(i=0;i<=N;i++){d=fabs(y[i]-x[i]);if(max<d) max=d;x[i]=y[i];}printf("%6d",k+1);for(i=1;i<=N;i++)printf("%12.8g",x[i]);printf("\n");k++;}while((max>=EPS)&&(k<NO));printf("\nk=%d\n",k);if(k>=NO) printf("\nfail!\n");elsefor(i=1;i<=N;i++)printf("x[%d]=%g\t",i,x[i]);}结果拉格朗日插值多项式#include<stdio.h>#include<math.h>#define N 4doublex[N]={0.56160,0.56280,0.56401,0.56521},y[N]={0.82741,0.82659,0.82577,0.82495}; void main(){double x=0.5635;double L(double xx);double lagBasis(int k,double xx);void output();output();printf("\n近似值L(%g)=%g\n",x,L(x));}double lagBasis(int k,double xx){double lb=1;int i;for(i=0;i<N;i++)if(i!=k) lb*=(xx-x[i])/(x[k]-x[i]);return lb;}double L(double xx){double s=0;int i;for(i=0;i<=N;i++)s+=lagBasis(i,xx)*y[i];return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<N;i++)printf("\t%g",x[i]);printf("\nyi:");for(i=0;i<N;i++)printf("\t%g",y[i]);}结果牛顿插值多项式#include <math.h>#include <stdio.h>int a;#define M 4double x[M+1]={0.4,0.55,0.65,0.8,0.9},y[M+1]={0.41075,0.57815,0.69675,0.88811,1.02652}; void main(){double x;printf("输入x=");scanf("%lf",&x);printf("次数：");scanf("%d",&a);double N(double xx,int a);void output();output();printf("\n%d次牛顿插值多项式N(%g)=%g\n",a,x,N(x,a));}double N(double xx,int a){double s=y[0],d=1;int i,j;double df[M+1][M+1];for(i=0;i<=M;i++)df[i][0]=y[i];for(j=1;j<=a;j++)for(i=j;i<=a;i++)df[i][j]=(df[i][j-1]-df[i-1][j-1])/(x[i]-x[i-j]);printf("\nx\tf(x)\t");for(j=1;j<=a;j++) printf("%5d阶",j);for(i=0;i<=a;i++){printf("\n%g\t%g",x[i],y[i]);for(j=1;j<=i;j++)printf("\t%7.5g",df[i][j]);}for(i=1;i<=a;i++){d*=(xx-x[i-1]);s+=df[i][i]*d;}return s;}void output(){int i;printf("\n各节点信息:\nxi:");for(i=0;i<=M;i++)printf("\t%7g",x[i]);printf("\nyi:");for(i=0;i<=M;i++)printf("\t%7g",y[i]);}结果复合梯形公式#include<stdio.h>#include<math.h>#define f(x) 1/(x*x+1)#define Pi 3.1415926void main(){double a=0,b=1;double T,h,x;int n,i;printf("please input n:");scanf("%d",&n);h=(b-a)/n;x=a;T=0;for(i=1;i<n;i++){x+=h;T+=f(x);}T=(f(a)+2*T+f(b))*h/2;printf("T(%d)=%g\n",n,T);printf("The exact value is %g\n",Pi/4);}复合辛普森公式#include<stdio.h>#include<math.h>#define f(x) 1/(1+x*x)#define Pi 3.1415926void main(){double a=0,b=1;double S,h,x;int n,i;printf("please input Even n:");scanf("%d",&n);h=(b-a)/n;x=a; S=0;for(i=1;i<n;i++){x+=h;if(i%2==0) S+=2*f(x);else S+=4*f(x);}S=(f(a)+S+f(b))*h/3;printf("S(%d)=%g\n",n,S);printf("The exact value is %g\n",Pi/4);}龙贝格公式加速#include<stdio.h>#include<math.h>#define f(x) sin(x)/(1+x)#define M 3void main(){double a=0,b=1;double Long(double a,double b);printf("近似值I=%g\n",Long(a,b));}double Long(double a,double b){int n=1,i=1,j=1;double T[M+1][M+1],h,x,sum;h=b-a;T[0][0]=(f(a)+f(b))/2;for(j=1;j<=3;j++){x=a;h/=2;n*=2;sum=0;for(i=1;i<=n;i+=2){x=a+i*h;sum+=f(x);}T[j][0]=T[j-1][0]/2+h*sum;}for(i=1;i<=M;i++)for(j=1;j<=i;j++){T[i][j]=(pow(4,j)*T[i][j-1]-T[i-1][j-1])/(pow(4,j)-1);}printf("k\tT0\tT1\tT2\tT3\n");for(i=0;i<=M;i++){printf("%d",i);for(j=0;j<=i;j++)printf(" %g",T[i][j]);printf("\n");}return T[M][M];}。

计算方法实验报告

2019年计算方法（B）实验报告姓名：学号：专业：课程：计算方法（B）目录一、实验综述 (1)二、实验内容 (1)2.1 实验一 (1)2.2 实验二 (2)2.3 实验三 (3)2.4 实验四 (4)2.5 实验五 (6)三、思考总结 (7)附件A1 (8)附件A2 (9)附件A3 (10)附件A4 (12)附件A5 (14)一、实验综述计算方法在工程实践中得到了广泛的应用，是理工类研究生必备的知识技能。

按照2019年计算方法课程学习要求，本文对计算方法上机题目进行了算法设计、分析，利用matlab 2019b版本对算法进行实现，最终形成了实验报告。

以下为本次实验报告具体内容，包括五个实验部分和一个思考总结部分。

二、实验内容2.1 实验一2.1.1 实验题目用Jacobi迭代和Gauss-Seidel迭代解电流方程组，使各部分电流的误差均小于10-3。

2.1.2 算法分析a)首先列出方程组的系数矩阵A以及等式右端的矩阵b，A=[28,-3,0,0,0;-3,38,-10,0,-5;0,-10,25,-15,0;0,0,-15,45,0;0,-5,0,0,30 ];b=[10;0;0;0;0];为了验证A是否收敛，我们通过判断系数矩阵A是否为严格对角占优矩阵进行确定。

如果是，则可以进行Jacobi迭代和Gauss-Seidel迭代（利用matlab程序验证后，证明了矩阵A为严格对角占优矩阵）；如果不是，则需要采用其他方法进行判断迭代是否收敛。

b)对矩阵A分裂成三部分，，其中D为A的对角矩阵，E为A的下三角矩阵的相反数，F为A的上三角矩阵的相反数。

c) Jacobi迭代。

取x得初始向量为x=[0;0;0;0;0]，利用迭代公式进行循环计算，当的无穷范数小于10-3，即，停止循环。

d) Gauss-Seidel迭代。

取x得初始向量为x=[0;0;0;0;0]，利用迭代公式进行循环计算，当的无穷范数小于10-3，即，停止循环。

计算方法作业

应用计算方法学院：自动化学院班级：硕1104班姓名：陈兆辉计算方法2班教师：张晓丹目录上机作业1 (3)上机作业2 (6)上机作业3 (10)上机作业4 (12)上机作业5 (17)上机作业6 (18)上机作业7 (19)上机作业11.分别用不动带你迭代法与Newton 法求解方程2x- ex +3=0的正根和负根解：（1）用不动点迭代法求方程的正根1）迭代公式：()231ln k P k P ++=，初值p0=1.02）程序：k=1;Tol=0.00001;p0=1.0;N=100;while k<=Np=log(2*p0+3);if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3)显示结果：P= 1.9239K=11（2）用不动点迭代法求解方程的负根1）迭代公式： 1(3)/2k p k p e +=-，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k<=Np=(exp(p0)-3)/2;if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P= -1.3734k=7(3)用牛顿法求方程的正根1）迭代公式：()1232k k p k k k p p e p p e +-+=--，初值p0=1.02）程序：k=1;Tol=0.00001;p0=1.0;N=100;while k<=Np=p0-(2*p0-exp(p0)+3)/(2-exp(p0));if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P=1.9239K=8（4）用牛顿法求解方程的负根1）迭代公式：()1232k k p k k k p p e p p e +-+=--，初值p0=-1.02）程序：k=1;Tol=0.00001;p0=-1.0;N=100;while k<=Np=p0-(2*p0-exp(p0)+3)/(2-exp(p0));if abs(p-p0)<Tolbreak;endk=k+1;p0=p;enddisp(p);disp(k)3）显示结果：P= -1.3734K= 44）结果分析：不动点迭代法是求解非线性方程的主要方法，其中牛顿法应用最广泛，它求解方程的单根时具有二阶收敛性，但是它要求较好的初始近似值，牛顿法比普通的迭代法收敛速度快，能较少的迭代达到理想的值。

计算方法上机作业

西安交通大学计算方法B大作业

西安电子科技大学出版社计算方法上机答案

东南大学计算方法与实习上机实验一

应用计算方法

计算方法与计算 实验一误差分析

计算方法大作业1 克服Runge现象

计算方法上机程序

( 鞋盒)上机纸张大小计算方法

《计算方法》上机实验指导书刘轶中-8页word资料

计算方法上机作业

数值计算方法上机实验报告

数值计算方法上机题

数值计算方法上机实验题

计算方法大作业——三次样条插值

数值计算方法上机实习题答案.doc

中国石油大学(华东)计算方法上机实验程序

计算方法实验上机报告(完整版)

计算方法实验报告

计算方法作业

计算方法与计算实验一误差分析