第10章工程数值计算方法实验指导

1、第10章 工程数值计算方法实验指导10.1 引言数值计算主要用于解决工程中常见的基本数学问题，研究其相关的数值解法，包含了线性方程组、非线性方程组的解法，矩阵求逆、数值逼近、数值微积分、常微分方程及偏微分方程的数值解法等。数值计算的基本理论和研究方法建立在数学建模的基础上，与计算机科学密切相关。采用计算机解决这些问题的主要步骤是：分析实际问题、建立数学模型、设计算法、编写程序代码并上机实现。本章精选前面章节中重点内容组成了7个实验项目，要求学生上机实现各个实验项目并编写实验报告。实验一 线性方程组的直接解列主元消去法解线性方程组；实验二 线性方程组的迭代解雅可比法、高斯赛德尔迭代法解线性方程组

2、；实验三 非线性方程的近似解二分法、牛顿法求非线性方程的根；实验四 插值问题拉格朗日插值、牛顿插值；实验五 曲线拟合问题最小二乘法；实验六 数值积分复化辛甫生公式；实验七 求解常微分方程的初值问题改进欧拉方法、四阶龙格库塔方法。10.2 实验一 线性方程组的直接解列主元消去法解线性方程组一、实验目的掌握列主元消去法解线性方程组的理论。二、实验环境PC机一台，C语言、MATLAB任选。三、实验内容用Gauss列主元消去法求解方程组。四、实验原理Gauss列主元消去法第一步：消元过程对于,将进行变换为，其中是上三角矩阵。即(1) 选列主元 选取第列中绝对值最大元素作为主元。(2) 换行(3) 归一

3、化(4) 消元第二步：回代过程由解出。五、实验步骤1要求上机实验前先编写出程序代码； 2编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、思考题 进一步思考全主元消去法的实现步骤。七、参考代码#include stdio.h #include math.h #define n 3 main() int i,j,k; int mi; float mv,tmp; float ann=0.01,2,-0.5,-1,-0.5,2,5,-4,0.5; float bn=-5,5,9,xn; for(k=0;kn-1;k+) mi=k; mv

4、=fabs(akk); for(i=k+1;imv) mi=i; mv=fabs(aik); if(mik) tmp=bk; bk=bmi; bmi=tmp; for(j=k;jn;j+) tmp=akj; akj=amij; amij=tmp; for(i=k+1;in;i+) tmp=aik/akk; bi=bi-bk*tmp; for(j=k+1;j=0;i-) xi=bi; for(j=i+1;jn;j+) xi=xi-aij*xj; xi=xi/aii; printf(nThe result is:); for(i=0;in;i+) printf(nx%d=%4.2f,i,xi);

5、程序运行后，输出结果如下10.3 实验二 线性方程组的迭代解雅可比迭代法、高斯赛德尔迭代法解线性方程组一、实验目的掌握雅可比迭代法和高斯赛德尔迭代法求解线性方程组的理论。二、实验环境PC机一台，C语言、MATLAB任选。三、实验内容分别采用雅可比迭代法、高斯赛德尔迭代法求解线性方程组。四、实验原理1雅可比迭代法设方程组的系数矩阵的对角线元素，为迭代次数容许的最大值，为容许误差。(1) 取初始向量，令。(2) 对，计算(3) 如果 并且，则输出方程的解，结束；如果，则说明方程组不收敛，输出方程组无解，终止程序；否则，转(2)。2高斯-赛德尔迭代法(1) 判断线性方程组是否主对角占优，即判断是否满

6、足。(2) 取初始向量，令。(3) 直接分离，并建立高斯-赛德尔迭代格式五、实验步骤1要求上机实验前先编写出程序代码； 2编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、实验注意事项雅可比迭代法和高斯赛德尔迭代法解线性方程组应注意迭代初值的选取和迭代终止的条件。七、思考题为节省计算时间，提高精度，如何选择合适的数据结构使得迭代收敛的速度更快？八、参考代码1雅可比迭代法#include#include#include#include#define EPS 1e-6#define MAX 100float * Jacobi(floa

7、t a34,int n)float *x,*y,epsilon,s;int i,j,k=0;x=(float *)malloc(n*sizeof(float);y=(float *)malloc(n*sizeof(float);for ( i = 0; i n; i+)xi=0;while(1)epsilon=0;k+;for (i = 0; i n; i+)s=0;for (j = 0; j n; j+)if(j=i)continue;s+=aij*xj;yi=(ain-s)/aii;epsilon+=fabs(yi-xi);if (epsilon=MAX) printf(The Metho

8、d is disconvergent!n);return y;for (i = 0; i n; i+)xi=yi; main()int i;float a34=5,2,1,8,2,8,-3,21,1,-3,-6,1;float *x;x=(float *)malloc(3*sizeof(float);x=Jacobi(a,3);for ( i = 0; i 3; i+)printf(x%d=%fn,i,xi );getch();程序运行后，输出结果如下2高斯赛德尔迭代法#include stdio.h #include math.h #define MAX 100 #define n 3 #d

9、efine exp 0.005 main() int i,j,k,m; float temp,s; float ann=5, 2,1,2,8,-3,1,-3,-6; float static bn=8,21,1; float static xn=0,0,0,Bnn,gn; for(i=0;in;i+) for(j=0;jn;j+) Bij=-aij/aii; gi=bi/aii; for(i=0;in;i+) Bii=0; m=0; do temp=0; for(i=0;in;i+) s=xi; xi=gi; for(j=0;jtemp) temp=fabs(xi-s); m+; printf

10、(n%dth result is:,m); printf(nx0=%7.5f,x1=%7.5f,x2=%7.5f,x0,x1,x2); printf(ntemp=%f,temp); while(temp=exp); printf(nnThe last result is:); for(i=0;in;i+) printf(nx%d=%7.5f,i,xi); 输出结果如下10.4 实验三 非线性方程的近似解二分法、牛顿法求非线性方程的根一、实验目的掌握二分法与牛顿迭代法求非线性方程组的根的理论。二、实验环境PC机一台，C语言、MATLAB任选。三、实验内容用二分法计算在上的单根，要求精确到小数点后

11、4位。用牛顿迭代法计算的近似解，要求其误差不超过。四、实验原理二分法通过将含根区间逐步二分，从而将根的区间缩小到容许误差范围内；牛顿法通过迭代逐步求解出近似解。 1二分法设为连续函数，方程的隔根区间为，设，。首先将区间二分得中点，计算在中点的函数值。若，则就是方程的根；否则，若，则隔根区间变为。若，则方程的有根区间变为，将新的隔根区间记为。再将二分，重复上述过程，就得到一系列隔根区间并有，。当时，趋近于零，收敛于点，此点即为方程的根。2牛顿法设是方程隔根区间中的一个初始值，把在点附近展开成泰勒级数取其线性部分作为非线性方程的近似方程，则有设，则其解为再把在附近展开成泰勒级数，也取其线性部分作为

12、 的近似方程。若，则得由此得到牛顿迭代法的一个迭代序列 五、实验步骤1要求上机实验前先编写出程序代码； 2编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、实验注意事项二分法应注意区间二分后，判断根可能的区间范围，牛顿法应注意确定迭代过程的结束的判定。七、参考代码1二分法#include#include #define f(x) (x*x-1)*x-1) void main() float a,b,x,eps; int k=0; printf(intput epsn); scanf(%f,&eps); printf(a,b=n);

13、 for(;) scanf(%f, %f,&a ,&b); if(f(a)*f(b)=0) printf(二分法不可使用,请重新输入:n); else break; do x=(a+b)/2; k+; if(f(a)*f(x)0) a=x; else break; while(fabs(b-a)eps); x=(a+b)/2; printf(n The root is x=%f, k=%dn,x,k);程序运行后，输出结果如下2牛顿法#include #include float f(float x) return(-3*x+4)*x-5)*x+6; float f1(float x) ret

14、urn (-9*x+8)*x-5; void main() float eps,x0,x1=1.0; printf(input eps:n); scanf(%f,&eps); do x0=x1; x1=x0-f(x0)/f1(x0); while(fabs(x1-x0)eps); printf(x=%fn,x1);程序运行后，输出结果如下10.5 实验四 插值问题拉格朗日插值、牛顿插值一、实验目的掌握拉格朗日插值法与牛顿插值法的理论。二、实验环境PC机一台，C语言、MATLAB任选。三、实验内容已知下列数据表0.561600.562800.564010.565210.827410.826590

15、.825770.82495分别采用拉格朗日插值公式和牛顿插值公式计算=0.5635时的函数值。四、实验原理已知个插值节点的函数值，分别采用由拉格朗日插值公式、牛顿插值公式构造插值多项式，并由该插值多项式计算待求点的函数值。1Lagrange插值公式输入，取次多项式作为插值函数，构造一组插值基函数它满足，并称为拉格朗日插值基函数。插值多项式为，其满足插值条件 ，并称此多项式为拉格朗日插值多项式。特别地，当时称为线性插值，其插值多项式为满足，从几何上看，为过两点的直线。当时，称为抛物线插值，其插值多项式为 满足。从几何上看为过点和的一条抛物线。2Newton插值公式输入，对，计算函数的各阶差商，计

16、算Newton插值函数值五、实验步骤1要求上机实验前先编写出程序代码； 2编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、思考题比较Lagrange插值法与Newton插值法的计算复杂性。七、参考代码1拉格朗日插值#include#include#includefloat Lagrange(float *x,float *y,float xx,int n)int i,j;float *a,yy=0.0;a=(float *)malloc(n*sizeof(float);for(i=0;i=n-1;i+)ai=yi;for(j=0

17、;j=n-1;j+)if(j!=i)ai*=(xx-xj)/(xi-xj);yy+=ai;free(a);return yy;void main()float x4=0.56160,0.56280,0.56401,0.56521;float y4=0.82741,0.82659,0.82577,0.82495;float xx=0.5635,yy;float Lagrange(float *,float *,float ,int);yy=Lagrange(x,y,xx,4);printf(x=%f,y=%fn,xx,yy);getch();程序运行后，输出结果如下2牛顿插值#include s

18、tdio.hmain() float static x4=0.56160,0.56280,0.56401,0.56521; float static y4=0.82741,0.82659,0.82577,0.82495; int i,k; double c,p; for (k=1;k=k;i-) yi=(yi-yi-1)/(xi-xi-k); printf(n%8.6f,yi); c=0.5635; printf(np%4.2f=%8.6f,c,y0+y1*(c-x0)+y2*(c-x0)*(c-x1)+y3*(c-x0)*(c-x1)*(c-x2)+y4*(c-x0)*(c-x1)*(c-x

19、2)*(c-x3); p=y4; for(i=3;i=0;i-) p=p*(c-xi)+yi; 程序运行后，输出结果如下10.6 实验五 曲线拟合问题最小二乘法一、实验目的掌握最小二乘法的基本原理。 二、实验环境PC机一台，C语言、MATLAB任选。三、实验内容已知一组数据00.51.21.62.22.633.64.203.24.24.86.16.98.29.912试用最小二乘法求出与此组数据相拟合的多项式函数。四、实验原理已知数据对，求多项式，使得为最小。注意到此时，多项式系数满足下面的线性方程组其中， ， 。调用解线性方程组的函数程序即可。五、实验步骤1要求上机实验前先编写出程序代码； 2

20、编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、参考代码用最小二乘法求拟合函数#include#include#define MAX 100void main() int i,j,k,m,n,N,mi; float tmp,mx; float XMAXMAX,YMAX,xMAX,yMAX,aMAX; printf(n 输入拟合多项式的次数:n); scanf(%d,&m); printf(n 输入给定点的个数n及坐标(x,y):n); scanf(%d,&N); printf(n); for(i=0;iN;i+) scanf(%

21、f,%f,&xi,&yi); for(i=0;i=m;i+) for(j=i;j=m;j+) tmp=0; for(k=0;kN;k+) tmp=tmp+pow(xk,(i+j); Xij=tmp; Xji=Xij; for(i=0;i=m;i+) tmp=0; for(k=0;kN;k+) tmp=tmp+yk*pow(xk,i); Yi=tmp; for(j=0;jm;j+) for(i=j+1,mi=j,mx=fabs(Xjj);imx) mi=i; mx=fabs(Xij); if(jmi) tmp=Yj;Yj=Ymi;Ymi=tmp; for(k=j;k=m;k+) tmp=Xjk;

22、 Xjk=Xmik; Xmik=tmp; for(i=j+1;i=m;i+) tmp=-Xij/Xjj; Yi+=Yj*tmp; for(k=j;k=0;i-) ai=Yi;for(j=i+1;j=m;j+) ai-=Xij*aj; ai/=Xii;printf(n 所求的二次多项式为:n);printf(P(x)=%f,a0);for(i=1;i=m;i+)printf(+(%f)*x%d,ai,i);printf(n);程序运行后，输出结果如下10.7 实验六 数值积分复化辛甫生公式一、实验目的掌握复化辛甫生公式计算定积分的理论。二、实验环境PC机一台，C语言、MATLAB任选。三、实验内

23、容用复化辛甫生公式求，要求误差不超过 。四、实验原理复化辛甫生公式为。五、实验步骤1要求上机实验前先编写出程序代码； 2编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、实验注意事项在积分中，注意构造函数下标的选取。七、参考代码采用复化辛甫生公式，计算定积分的值#include#includevoid main()int i,n=2;float s;float f(float);float Simpson(float (*)(float),float,float,int);for (int i = 0; i = 2; i+)s=S

24、impson(f,0,1,n);printf(s=%fn,s );getch();float Simpson(float (*f)(float),float a,float b,int n)int k;float s,s1,s2=0.0;float h=(b-a)/n;s1=f(a+h/2);for (int k = 1; k =n-1 ; k+)s1+=f(a+k*h+h/2);s2+=f(a+k*h);s=h/6*(f(a)+4*s1+2*s2+f(b);return s;float f(float x)return 4/(1+x*x);程序运行后，输出结果如下10.8 实验七 求解常微分

25、方程的初值问题改进欧拉方法、四阶龙格库塔方法一、实验目的掌握用改进欧拉法、四阶龙格-库塔法求解一阶常微分方程的初值问题的理论。二、实验环境PC机一台，C语言、MATLAB任选。三、实验内容1、用改进欧拉算法求微分方程 的数值解，步长。 2、用四阶龙格库塔法求微分方程的数值解，步长。四、实验原理1改进欧拉法公式2四阶龙格库塔法公式五、实验步骤1要求上机实验前先编写出程序代码； 2编辑录入程序；3调试程序并记录调试过程中出现的问题并修改程序；4记录运行时输入数据和输出结果；5撰写实验报告。六、参考代码1改进欧拉法计算微分方程#include #define MAXSIZE 20double f(double x,double y);void main(void)double a,b,h,xMAXSIZE,yMAXSIZE;long i,n;printf(n请输入求解区间a,b：);scanf(%lf,%lf,&a,&b);printf(n请输入步长h：);scanf(%lf,&h);n=(long)(b-a)/h);x0=a;printf(n请输入起点x0=%lf处的纵坐标y0：,x0);scanf(%