数值计算基础实验指导部分实验源代码复习指导三套试题及其答案.doc

数值计算基础实验指导书2010 年目录实验一 直接法解线性方程组的 . 3 实验二 插值方法. 12 实验三 数值积分. 6 实验四 常微分方程的数值解 . 8 实验五 迭代法解线性方程组与非线性方程 . 10实验一 直接法解线性方程组一、实验目的掌握全选主元消去法与高斯-塞德尔法解线性方程组。二、实验内容分别写出 Guass 列选主元消去法与追赶法的算法，编写程序上机调试出结果，要求所编 程序适用于任何一解线性方程组问题，即能解决这一类问题，而不是某一个问题。实验中以 下列数据验证程序的正确性。1、用 Guass 列选主元消去法求解方程组2.52.3- 5.1 x1 3.75.39.61.5 x = 3.88.11.7- 4.3 2 x3 5.52、用追赶法求解方程组- 20000 x1 - 10 1- 2000 x 0 2 01- 200 x3 = 0 001- 20 x4 0 0001- 2x50三、实验仪器设备与材料主流微型计算机四、实验原理1、Guass 列选主元消去法对于 AX =BA B 是上三角矩阵。即：1）、消元过程：将（A|B）进行变换为 ( | ) ，其中 A aaLab 1aLab 1112 1n1 12 1n 1 a21 a22 La2 nb2 01La2 nb2 MMMM MMMM an1an 2bLann n 00bLann n k 从 1 到 n-1 a、 列选主元选取第 k 列中绝对值最大元素 max aik k in作为主元。b、 换行akjbkc、 归一化 aij , j = k + 1,L, n bid、 消元akj / akkbk / akk akj , j = k + 1,L, n bkaij - aik akj aij , i = k + 1,L, n; j = k + 1,L, n bi - aik bk bi , i = k + 1,L, n2）、回代过程：由 ( | ) 解出 x , x,L, x 。A Bbn / ann xn nn n-1 1bk - akj x j xk , k = n - 1,L,2,1j =k +12、追赶法线性方程组为： ac x f 111 1 b2a2c2Ob3a3c3 x2 3 x f 2 3 fOOOO M = M bn -1Oan -1cn-1M xn-1M fn-1bnan xn f n 做 LU 分解为：a 1b 11 g 2a 2g a1b 2OOL = 33, R = OObOOO O1n-1 g na n 1 分解公式：g= a ii(i = 2,3,L, n)a1 = b1 ,a i = bi - g i b i -1(i = 2,3,L, n)ib = cia i(i = 1,2,L, n - 1)则Ax = f LUx = fLy = f Ux = y回代公式：y=f1 1y ia1= f i - g i yi -1a i(i = 2,3,L, n)xn = ynxi = yi - b i xi +1(i = n - 1, n - 2,L,1)五、实验步骤1、理解并掌握全选主元消去法与高斯-塞德尔迭代法公式；2、画出全选主元消去法与高斯-塞德尔迭代法的流程图3、使用 C 语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、 实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项注意如何定义数据结构以保存矩阵和解以降低算法的复杂性。八、思考题若使用全主元消去法，在编程中应如何记录保存对于未知数的调换。实验二 插值方法一、实验目的掌握拉格郎日插值法与牛顿插值法构造插值多项式。二、实验内容分别写出拉格郎日插值法与牛顿插值法的算法，编写程序上机调试出结果，要求所编程 序适用于任何一组插值节点，即能解决这一类问题，而不是某一个问题。实验中以下列数据 验证程序的正确性。已知下列函数表xi0.561600.562800.564010.56521yi0.827410.826590.825770.82495求 x=0.5635 时的函数值。三、实验仪器设备与材料主流微型计算机四、实验原理已知 n 个插值节点的函数值，则可由拉格郎日插值公式与牛顿插值公式构造出插值多项 式，从而由该插值多项式求出所要求点的函数值。拉格郎日插值公式与牛顿插值公式如下：1、Lagrange 插值公式nLn ( x) = l0 ( x) y0 + l1 ( x) y1 + . + ln ( x) yn = yk lk ( x)k =0( x - x)(x - x )L( x - x)(x - x)L( x - x )n x - xl ( x) = 0 1 k -1 k +1 n = jkk ( x -x0 )(xk -x1 )L( xk -xk -1)(xk -xk +1)L( xk -xn )j =0j kxk - x j2、Newton 插值公式N n ( x) =f ( x0 ) + f x0 , x1 (x - x0 ) + f x0 , x1 , x2 (x - x0 )(x - x1 )+ L + f x0 , x1 ,L xn (x - x0 )(x - x1 )L( x - xn-1 )五、实验步骤1、理解并掌握拉格郎日插值法与牛顿插值法的公式；2、画出拉格郎日插值法与牛顿插值法算法的流程图；3、使用 C 语言编写出相应的程序并调试验证通过。六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项Newton 插值法在编程时应注意定义何种数据结构以保存差商。八、思考题比较 Lagrange 插值法与 Newton 插值法的异同。实验三 数值积分一、实验目的掌握复化梯形法与龙贝格法计算定积分。二、实验内容分别写出变步长梯形法与 Romberge 法计算定积分的算法，编写程序上机调试出结果， 要求所编程序适用于任何类型的定积分，即能解决这一类问题，而不是某一个问题。实验中以下列数据验证程序的正确性。求 1sin x dx, e 0.00001 。0 x三、实验仪器设备与材料主流微型计算机四、实验原理通过变步长梯形法与龙贝格法，我们只要知道已知 n 个求积节点的函数值，则可由相应 的公式求出该函数的积分值，从而不需要求该函数的原函数。变步长梯形法与龙贝格法公式 如下：1、变步长梯形法 2n-1 hTn = f ( xi ) + f ( xi +1 )i =0hn-1= f (a) + 2 f ( xi ) + f (b)2i =11 h n-1T2 n = Tn + f ( xi +1 / 2 )2 2 i =0用 T2n - Tn e 来控制精度2、龙贝格法1h n-1T2 n = Tn + f ( xi +1 / 2 )2 2 i =0S n =4 1T2n - Tn3 3Cn =16 1S 2n - S n15 1564 1Rn =C2n - Cn63 63用 R2n - Rn e 来控制精度五、实验步骤1、理解并掌握变步长梯形法与龙贝格法的公式；2、画出变步长梯形法与龙贝格法的流程图3、使用 C 语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、 实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项1 sin x在 0 x的定义。dx 积分中，被积函数在 x=0 点函数值为 1，对该点在程序设计中应注意对其八、思考题使用复化梯形法与复化 Simpson 法来计算该问题有何缺点？实验四 常微分方程的数值解一、实验目的掌握改进欧拉法与四阶龙格-库塔求解一阶常微分方程的初值问题。二、实验内容分别写出改进欧拉法与四阶龙格-库塔求解的算法，编写程序上机调试出结果，要求所 编程序适用于任何一阶常微分方程的数值解问题，即能解决这一类问题，而不是某一个问题。 实验中以下列数据验证程序的正确性。 y = - xy 2求 y(0) = 2(0 x 5)步长 h=0.25。三、实验仪器设备与材料主流微型计算机四、实验原理常微分方程的数值解主要采用“步进式”，即求解过程顺着节点排列次序一步一步向前 推进，在单步法中改进欧拉法和四阶龙格-库塔法公式如下：1、改进欧拉法yn+1 = yn + hf (xn , yn )yn+1= yn+ h f ( x2 n, yn) + f ( xn+1, yn+1 )2、四阶龙格-库塔法h yn +1 = yn + 6 (k1 + 2k 2 + 2k3 + k 4 )k1 =k 2 =f ( xn , yn )f ( x + hn2h, yn +hk1 )2hk3 =k 4 =f ( xn +, yn +k 2 )22f ( xn + h, yn + hk3 )五、实验步骤1、理解并掌握改进欧拉法与四阶龙格-库塔法的公式；2、画出改进欧拉法与四阶龙格-库塔法的流程图3、使用 C 语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项 y = - xy 2 y(0) = 2精度的变化八、思考题(0 x 5)的精确解为 y = 2 /(1 + x 2 ) ，通过调整步长，观察结果的如何对四阶龙格-库塔法进行改进，以保证结果的精度。实验五 迭代法解线性方程组与非线性方程一、实验目的掌握高斯-塞德尔迭代法求解线性方程组与牛顿迭代法求方程根。二、实验内容分别写出高斯-塞德尔迭代法与牛顿迭代法的算法，编写程序上机调试出结果，要求所 编程序适用于任何一个方程的求根，即能解决这一类问题，而不是某一个问题。实验中以下 列数据验证程序的正确性。1、高斯-塞德尔迭代法求解线性方程组 721- 2 x1 4 9153- 2 x 7 2 = - 2- 2 115 x3 - 1 13213 x4 0 2、用牛顿迭代法求方程 x3 - x - 1 = 0 的近似根，e 0.00001，牛顿法的初始值为 1。三、实验仪器设备与材料主流微型计算机四、实验原理二分法通过将含根区间逐步二分，从而将根的区间缩小到容许误差范围。牛顿通过迭代 的方法逐步趋进于精确解，该两种方法的公式如下：1、高斯-塞德尔迭代法1）判断线性方程组是否主对角占优n aijj =1j i aii , i = 1,2,L, n2）直接分离 xi，即nxi = (di - bij x j ) / aii , i = 1,2,L, nj =1建立高斯-塞德尔迭代格式为：( k +1)xii -1= (d i - aij x j( k +1)n- aij x j( k ) / aii , i = 1,2,L, nj =1j =i +13）取初值迭代求解至所要求的精度为止。2、牛顿法kx = x- f ( xk )k +1k f ( x )五、实验步骤1、理解并掌握二分法与牛顿法的公式；2、画出二分法与牛顿法的流程图3、使用 C 语言编写出相应的程序并调试验证通过六、实验报告要求1、统一使用武汉科技大学实验报告本书写，实验报告的内容要求有：实验目的、 实验内容、程序流程图、源程序、运行结果及实验小结六个部分。2、源程序需打印后粘贴在实验报告册内；3、运行结果以屏幕截图形式保存并打印后粘贴在实验报告册内。七、实验注意事项对于二分法应注意二分后如何判断根的区间，对于牛顿法注意如何确定迭代过程的结 束八、思考题若使用牛顿法是发散的，如何对牛顿法进行改进以保证其收敛性。前三个实验的程序代码（C/C+）和运行结果截图#include #include #include using namespace std;Gauss 全选主元解方程组的源程序及运行结果class Matrixpublic: Matrix();Matrix();void SetMatrix(const int n,const double esp1);/构造线性方程组相应的矩阵，n为方程的未知数数目，esp1 为要求的精度void Max(const int r);/全选主元void ChangeRC(const int r);/根据主元变换矩阵的行或列void Eliminate(const int r);/处理消元工作void Result()const;/计算方程的解void Calculate();int GetRank()const;/返回矩阵的行数double GetX(const i)const;/确定方程组的第 i 个解(1=i=N)private:/指针 a 和 b 分别用于存储方程组的未知数系数和方程=右边的常数，esp 存/储精度double *a,*b,esp;/指针 flag 用于记录方程组解的顺序int *flag;/以下的结构体用于在全选主元中记录最大主元的位置struct coordinateint row,column;location;int N;/方程组未知数的数目;int main()int n;double esp1; Matrix matrix; docoutn;if(n0) n=0;/如何控制非法字符的输入？while(n=0);docoutesp1;if(esp10) esp1=0;/输入不合法的精度就把精度置 0while(esp1=0);cout输入线性方程组的增广矩阵:n; matrix.SetMatrix(n,esp1);/设置矩阵内的数据 matrix.Calculate();/计算方程组的解/输出方程组的解coutnn 方程组的解如下:n;for(int i=1;i=matrix.GetRank();+i)coutXi:matrix.GetX(i)endl;return 0;Matrix:Matrix()/将 Matrix 类的数据成员初始化a=NULL; b=NULL; flag=NULL; location.row=0; location.column=0; esp=0;N=0;Matrix:Matrix()/释放指针 a、b 和 flag 指向的内存空间deletea;deleteb;deleteflag;void Matrix:SetMatrix(const int n,const double esp1)N=n;esp=esp1;a=new doubleN*N; b=new doubleN; flag=new intN;/判断是否成功分配存储区if(a=NULL|b=NULL|flag=NULL)cout分配存储区失败！n;exit(EXIT_FAILURE);/读取线性方程组的增广矩阵for(int i=0;iN;+i)for(int j=0;j*(a+i*N+j);cin*(b+i);/flag 中存储的值对应相应的 x 值，当方程的解由于列变换交换后，flag 中/的值也相应交换，最后用于恢复解的顺序for( i=0;iN;+i)*(flag+i)=i;void Matrix:Max(const int r)double max=0;for(int i=r;iN;+i)for(int j=r;jN;+j)if(maxfabs(*(a+i*N+j)max=fabs(*(a+i*N+j);/设定最大主元的行、列location.row=i;location.column=j;/最大主元小于输入的精度时，认为方程组无解，退出程序if(max=esp)cout方程组无解！n;exit(EXIT_FAILURE);void Matrix:ChangeRC(const int r)double temp;/如果最大主元所在的行不在当前行，则进行行变换if(location.row!=r)for(int i=r;iN;+i)temp=*(a+r*N+i);*(a+r*N+i)=*(a+location.row*N+i);*(a+location.row*N+i)=temp;temp=br; br=blocation.row; blocation.row=temp;/若最大主元所在的列不在当前的 r 列，则进行列变换if(location.column!=r)for(int i=r;iN;+i)temp=*(a+i*N+r);*(a+i*N+r)=*(a+i*N+location.column);*(a+i*N+location.column)=temp;/交换 flag 中的元素来标记方程解的位置变化int temp1;temp1=*(flag+r);*(flag+r)=*(flag+location.column);*(flag+location.column)=temp1;void Matrix:Eliminate(const int r)if(fabs(*(a+N*r+r)=esp)cout方程组无解！n;exit(EXIT_FAILURE);for(int i=r+1;iN;+i)for(int j=r+1;jN;+j)(*(a+i*N+j)-=(*(a+i*N+r)*(*(a+r*N+j)/(*(a+r*N+r); (*(b+i)-=(*(b+r)*(*(a+i*N+r)/(*(a+r*N+r);void Matrix:Result()constif(fabs(*(a+N*(N-1)+N-1)=esp)cout=0;-i)temp=0;for(int j=i+1;jN;+j)temp+=(*(a+i*N+j)*(*(b+j);*(b+i)=(*(b+i)-temp)/(*(a+i*N+i);/根据 flag 中的数据用冒泡排序法恢复方程组解的次序for(i=0;iN-1;+i)for(int j=0;j*(flag+j+1)int temp1;/交换解的顺序temp=*(b+j);*(b+j)=*(b+j+1);*(b+j+1)=temp;/交换用于标记的元素的顺序temp1=*(flag+j);*(flag+j)=*(flag+j+1);*(flag+j+1)=temp1;void Matrix:Calculate()/根据矩阵行数重复进行寻找最大主元、变换行或列、消元for(int i=0;iGetRank()-1;+i) Max(i); ChangeRC(i); Eliminate(i);Result();int Matrix:GetRank()constreturn N;/返回矩阵的行数double Matrix:GetX(const int i)constreturn *(b+i-1);运行结果追赶法求解方程组的算法：1 输入方程组的维数 n，将主对角元素 b(i)(i=0:n-1)，主对角元素左边的元素a(i)(i=0:n-2)，主对角元素右边的元素 c(i)(i=0:n-2)，右端项的元素 f(i)(i=0:n-1)2 对方程组的系数矩阵作 Crout 分解, (0)=b(0),对于 i=0:n-2,c(i):= (i):= c(i)/b(i), b(i+1):= (i+1):= b(i+1)-a(i)* (i)3.解方程组 Ly=fb(0):=y(0):=f(0)/ (0):=f(0)/b(0)对于 i=1：n-1，b(i):=y(i):=f(i)-a(i-1)*y(i-1)/b(i)4.解方程组 Ux=ya(n-1):=x(n-1):=b(n-1)对于 i=n-2:0,a(i)=x(i):=b(i)-c(i)*a(i+1);5输出方程组的解 a(0:n-1)#include#include用追赶法求解方程组的源程序及运行结果using namespace std;class MatrixThrpublic: MatrixThr();MatrixThr();void SetMatrixThr(const int n);/设置三对角矩阵的数据void Result();/计算三对角矩阵的解double GetX(const int i)const;/取得第 i 个解,i 从 1 开始int GetN() const;/返回未知数的数目private:int N;/N 为未知数的数目/b 为矩阵主对角线的元素首地址，a 为主对角线左边一斜条元素的首地址，/c 为主对角线右边一斜条元素首地址，f 为方程组的常数首地址double *a,*b,*c,*f;int main()MatrixThr matrix;int n;docoutn;while(n3);cout请依次输入三对角方程组每行的数据(0 元素除外):n; matrix.SetMatrixThr(n); /计算方程组的解 matrix.Result();/输出方程组的解cout方程的解如下:n;for(int i=1;i=matrix.GetN();+i)coutXi:matrix.GetX(i)endl;return 0;MatrixThr:MatrixThr()/初始化相关数据N=0; a=NULL; b=NULL; c=NULL;f=NULL;MatrixThr:MatrixThr()/释放分配的内存空间delete a; delete b; delete c;delete f;void MatrixThr:SetMatrixThr(const int n)/根据输入的未知数个数设置矩阵的数据N=n;a=new doubleN; b=new doubleN; c=new doubleN; f=new doubleN;/若内存分配失败，退出程序if(a=NULL|b=NULL|c=NULL|f=NULL)cout*b*c*f;for(int i=1;i*(a+i-1)*(b+i)*(c+i)*(f+i);cin*(a+i-1)*(b+i)*(f+i);void MatrixThr:Result()/对系数矩阵 A 作 Crout 分解for(int i=0;iN-1;+i)/将 U 中的 存于指针 b 中，L 中的 存于指针 c 中*(c+i)/=(*(b+i);*(b+i+1)-=(*(a+i)*(*(c+i);/解方程组 Ly=f，求得的 y 值存于指针 b 中*b=(*f)/(*b);for(i=1;i=0;-i)*(a+i)=*(b+i)-(*(c+i)*(*(a+i+1);int MatrixThr:GetN()constreturn N;double MatrixThr:GetX(const int i)constreturn *(a+i-1);运行结果#includeLagrange 插值法的源程序：#includeusing namespace std;class Lagrangepublic: Lagrange();Lagrange();void SetLagrange(const int n);/根据用户的输入设置 Lagrange 类中的插值点数据bool Exist(const double x,const int i);/检测是否输入了与前 i 个插值结点横坐标相同的点intGetN()const;/获取插值结点的数目void Calculate(const double a);/计算横坐标 a 对应的函数值double GetResult()const;/返回计算的函数值private:int N;/插值结点的数目double *x,*y,zx,zy;/x、y 分别用于存储插值点的数据，zx、zy 表示所求的坐标点;int main()int n=2; Lagrange L; doif(n2) cout2endl;coutn;while(n2);cout请输入各插值点横、纵坐标的数据n; L.SetLagrange(n);double a;couta; L.Calculate(a);cout横坐标a-对应的函数值为:L.GetResult()endl;return 0;Lagrange:Lagrange()/初始化相关数据N=0;x=y=NULL;zx=zy=0;Lagrange:Lagrange()/释放分配给指针的内存空间delete x;delete y;bool Lagrange:Exist(const double a,const int i)/遍历以输入的插值点，看是否重复插入横坐标相同的插值点for(int j=0;ji;+j)if(a=*(x+j) return true;return false;void Lagrange:SetLagrange(const int n)N=n;x=new doubleN;y=new doubleN;/判断是否成功为指针分配了内存空间if(x=NULL|y=NULL)cout分配存储空间失败！endl;exit(EXIT_FAILURE);/输入插值点for(int i=0;i*(x+i)*(y+i);/如果不是输入第一个坐标值，则会对输入的横坐标进行合法性检测while(i!=0&Exist(*(x+i),i)=true)cout输入了重复的横坐标-请重新输入第i+1*(x+i)*(y+i);void Lagrange:Calculate(const double a)zx=a;for(int i=0;iGetN();+i)/计算插值基函数double temp=1;for(int j=0;ji;+j) temp*=(zx-*(x+j)/(*(x+i)-*(x+j);for(+j;jGetN();+j) temp*=(zx-*(x+j)/(*(x+i)-*(x+j);zy+=(*(y+i)*temp);int Lagrange:GetN()constreturn N;/返回插值点的数目double Lagrange:GetResult() constreturn zy;/返回求得的函数值Lagrange 插值法的运行结果#include using namespace std;Newton 插值法的源程序：class Matrixpublic: Matrix();Matrix();void SetMatrix(const int n);/根据用户输入的插值点的数据设置计算结果的矩阵int GetN()const;/返回插值点的数目void Calculate();/计算差商double GetResult(double x)const;/根据输入的横坐标求函数值，返回运算结果private:double *a,*f;/a 和 f 分别用于存储插值点的横坐标和相应的函数值int N;/记录插值点的数目;int main()Matrix matrix;int n;docoutn;while(n2);cout输入插值点的数据：endl;matrix.SetMatrix(n);matrix.Calculate();double x;coutx;cout所求的函数值为：matrix.GetResult(x)endl;return 1;Matrix:Matrix()/初始化数据a=f=NULL;N=0;Matrix:Matrix()/释放指针 a、f 指向的内存空间delete a;delete f;void Matrix:SetMatrix(int n)N=n;/为插值点创建动态数组a=new doubleN;f=new doubleN;/输入插值点的数据for(int i=0;i*(a+i)*(f+i);int Matrix:GetN() constreturn N;void Matrix:Calculate()/将差商存储在一个一维数组内for(int i=0;ii;-j)*(f+j)=(*(f+j)-*(f+j-1)/(*(a+j)-*(a+j-i-1);double Matrix:GetResult(double x) const/利用差商和插值点的横坐标及第一个插值点的纵坐标计算函数值double result=*f;/指针 f 指向第一个插值点的纵坐标for(int i=1;iGetN();+i)double temp=1;for(int j=0;ji;+j)temp*=(x-*(a+j);/计算(x-x0)*(x-x1)*(x-x(i-1)result+=(*(f+i)*temp;/指针(f+i)指向第 i 号差商return result;Newton 插值法的运行结果变步长梯形法求定积分的源程序：#include#includeusing namespace std;double function(const double x);/求被积函数的值并返回/accumulate()为求定积分的函数，a、b 分别为积分的上下限,默认精度为 0.00001 double accumulate(const double a,const double b,const double eps=0.00001);int main()double a,b,eps;/a,b 分别为定积分的上限和下限,h 为步长,eps 为要求的精度a=0; b=1; eps=0.00001;cout(sinX/X)在 0 到 1 上的积分为:accumulate(a,b)eps)n*=2; h/=2;/步长折半 T1=T2;/利用