计算方法上机作业

数值计算方法I 上机实验考试题（两题任选一题）
１．小型火箭初始质量为900千克，其中包括600千克燃料。

火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。

设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2.
A. 建立火箭升空过程的数学模型（微分方程）；
B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度．
２．小型火箭初始质量为1200千克，其中包括900千克燃料。

火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。

设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为
0)0(,0,01222==≤≤-+⎪⎭
⎫ ⎝⎛-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻．
今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））：
现有两种估计比例系数k 的方法：
1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。

2．用这组数据拟合一个k ．
请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）．。

计算方法上机作业——龙格库塔法龙格库塔法（Runge-Kutta method）是一种常用于求解常微分方程（Ordinary Differential Equation，ODE）的数值解法。

它是由德国数学家卡尔·龙格（Carl Runge）和马丁·威尔海姆·库塔（Martin Wilhelm Kutta）在20世纪初提出的。

龙格库塔法的基本思想是通过数值逼近来计算微分方程的近似解。

在讲解龙格库塔法之前，我们先来简单回顾一下ODE的一阶常微分方程的基本形式：y′(y)=y(y,y)其中，y(y,y)是已知函数。

龙格库塔法的核心是使用差分逼近计算函数的斜率。

假设我们要求解的方程为：y′(y)=y(y,y),y(y)=y₀所需计算的点为y₀,y₁,...,yy，对应的函数值为y₀,y₁,...,yy，其中y是步长的个数。

龙格库塔法通过递推关系式来计算估计值，并不断更新当前点的函数值。

接下来以龙格库塔法的经典四阶形式为例进行说明。

该方法的基本方程如下：yy+1=yy+(y₁+2y₂+2y₃+y₄)/6y₁=ℎy(yy,yy)y₂=ℎy(yy+ℎ/2,yy+y₁/2)y₃=ℎy(yy+ℎ/2,yy+y₂/2)y₄=ℎy(yy+ℎ,yy+y₃)其中y表示当前步骤，ℎ表示步长，yy表示当前点的函数值，y₁,y₂,y₃和y₄则表示对应的斜率。

使用龙格库塔法，我们可以通过不断递归计算来求得指定区间（例如[y,y]）上的函数值。

具体步骤如下：1.确定求解区间[y,y]和初始点(y₀,y₀)以及步长ℎ。

2.初始化：设置yy=y₀，yy=y₀。

3.对所有y=0,...,y−1：计算y₁,y₂,y₃和y₄，根据上述递推关系式。

根据递推关系式计算yy+1更新当前点的函数值，即yy+1=y(yy+1)。

更新当前点的y值，即yy+1=yy+ℎ。

4.返回结果：最终求得的函数值。

需要注意的是，选择适当的步长对最终结果的精度和计算效率都有重要影响。

计算方法上机报告姓名：学号：班级：目录题目一------------------------------------------------------------------------------------------ - 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个有效数字，则又将如何进行计算。

数值计算方法上机实验报告
一、实验目的
本次实验的主要目的是熟悉和掌握数值计算方法，学习梯度下降法的
原理和实际应用，熟悉Python语言的编程基础知识，掌握Python语言的
基本语法。

二、设计思路
本次实验主要使用的python语言，利用python下的numpy，matplotlib这两个工具，来实现数值计算和可视化的任务。

1. 首先了解numpy的基本使用方法，学习numpy的矩阵操作，以及numpy提供的常见算法，如矩阵分解、特征值分解等。

2. 在了解numpy的基本操作后，可以学习matplotlib库中的可视化
技术，掌握如何将生成的数据以图表的形式展示出来。

3. 接下来就是要学习梯度下降法，首先了解梯度下降法的主要原理，以及具体的实际应用，用python实现梯度下降法给出的算法框架，最终
可以达到所期望的优化结果。

三、实验步骤
1. 熟悉Python语言的基本语法。

首先是熟悉Python语言的基本语法，学习如何使用Python实现变量
定义，控制语句，函数定义，类使用，以及面向对象编程的基本概念。

2. 学习numpy库的使用方法。

其次是学习numpy库的使用方法，学习如何使用numpy库构建矩阵，学习numpy库的向量，矩阵操作，以及numpy库提供的常见算法，如矩阵分解，特征值分解等。

3. 学习matplotlib库的使用方法。

西安电子科技大学出版社《计算方法》任传祥等编著第九章计算方法上机参考答案实验一，算法一#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++程序进行计算验证采用从大到小或者从小到大的递推公式算法得到的数值基本稳定且误差不大。

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。

《线性代数》上机作业（一）
课 程《线性代
数》
上 机
内 容
向量与矩阵运算成 绩
姓 名王伟专 业
班 级
资源勘查工程2010级2
班
学 号1002040316
教学班指 导
教 师
熊文军
上 机
日 期
2011-10-20
一、上机目的
1、突出教学的重点，减少教学的难点；
2、掌握常用计算方法和处理问题的方法.
二、上机内容
1、向量与矩阵运算。

3、 上机作业
1、试分别生成 4 阶的单位矩阵和5 阶均匀分布的随机矩阵；
4 阶的单位矩阵:
5 阶均匀分布的随机矩阵:
2、生成列向量x=[1, 3, 5, 7, 9, … , 15]；
3、设A=[2 1 2 4; 1 2 0 2; 4 5 2 0; 0 1 1 7];
B=[12 1 2 4; 11 2 0 2; 14 5 2 0; 10 1 1 7];
求A+B，A-B，A*B，A\B，5*B，A^5，A.*B， A.^5及矩阵A的行列式，特征值，特征向量，秩和行最简形.
A+B
A-B
A*B
A\B
5*B
A^5
A.*B
A.^5
A的行列式
A的秩
A的行最简形
A的特征向量v和特征值 d
4、 上机心得体会。