贵州大学-数值分析上机

1、数值分析上机实验报告 数值分析上机实验报告学 院： 土木工程学院 专业名称： 土木工程 姓 名： XXXX 学 号： 201XXXXXXX 指导老师： XXX老师 2017年01月02日20第一题 一，用Newton法求方程：x7-28x4+14=0 在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001.)解：（1）解题方法及理论依据：Newton迭代法定理：设函数在有限区间上二阶导数存在，且满足条件： ；在区间上不变号；,其中是中使达到的一个。则对任意初始近似值，由Newton迭代过程 所生的迭代序列平方收敛于方程在区间上的唯一解。(2) 计算机程序C语言

2、#include#includevoid main()int flag=1,c; /设置flag控制while循环（即控制整个程序是否继续运行）while(flag=1) float x,y,f,f1; int i=0; printf(n请输入x的值：n); scanf(%f,&x); /输入x的初值 for(;(f1!=0)&(fabs(y-x)=0.000001);) /迭代循环控制条件 y=x; f=pow(x,7)-28*pow(x,4)+14; /原方程 f1=7*pow(x,6)-28*4*pow(x,3); /原方程的导数方程 x=x-f/f1; /牛顿迭代公式 i+; /累计迭

3、代的次数 printf(x=%f,y=%fn,x,y); /输出每次迭代的值 printf(n方程的近似根:y=%fn,y); printf(n迭代次数:%dn,i); printf(是否继续(Y/N)?); /提示是否输入下一个x的初值 getchar(); c=getchar(); if(c=N|c=n) /如果输入N或n，则flag=0退出程序；则继续运行程序 flag=0; else flag=1; 实验结果（截屏图）：（3）Matlab语言xk=1.2;fx=xk.7-28*xk.4+14;fxx=7*xk.6-4*28*xk.3;iffxx=0.0dispendxkk=xk-fx/

4、fxx;while abs(xkk-xk)=0.000001xk=xkk;fxx=7*xk.6-4*28*xk.3;fx=xk.7-28*xk.4+14;iffxx=0.0dispendxkk=xk-fx/fxx;endxkk实验结果（截屏图）：问题讨论与分析：Newton迭代法是一个二阶收敛迭代式，他的几何意义是的切线与x轴的交点,故也称为切线法。它是平方收敛的，但它是局部收敛的，即要求初始值与方程的根充分接近，所以在计算过程中需要先确定初始值。本题在讨论中，讨论了区间(0.1,1.9)两端点是否能作为Newton迭代的初值，结果发现0.1不满足条件，而1.9满足，能作为初值。第二题二，已知

5、函数值如下表：x1 2345f(x)00.6931471.09861231.38629441.6094378x678910f(x)1.79175951.94591012.0794452.19722462.3025851f(x)f(1)=1f(10)=0.1试用三次样条插值求f(4.563)和f(4.563)的近似值。解：（1）解题方法及理论依据： 公式1 公式2三弯矩方程：其矩阵形式：其中对于上述矩阵方程，由于其系数矩阵为三对角矩阵，所以用追赶法来解。解三对角方程组的追赶法：设线性方程组Ax=b的系数矩阵A为三对角矩阵,这时A可以作如下LU分解：(1) LU分解：首先，对i=2，3，.,n,计

6、算 ,(2) 解Lz=b：首先,对i=2，3，.,n,计算 (3) 解Ux=z：首先 ,对i=n-1,n-2，.,2,1,计算 最后把追赶法求得的唯一解代入上面的公式1和公式2，即可求得。（2） 计算机程序C语言#include #include void main()float X10=1,2,3,4,5,6,7,8,9,10; /输入节点float f10=0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.94594101,2.079445,2.1972246,2.3025851 ; /输入节点处的函数值float h10; fl

7、oat u9,v9,c9;float r10,z10,l10,a10,d10,m10;int i;float x=4.563;float S,s;for(i=0;i=8;i+)hi=Xi+1-Xi; /计算步长hd0=6*(f0-f1)/(X0-X1)-1)/h1; d8=6*(0.1-(f8-f9)/(X8-X9)/h8;u8=1;v0=1;for(i=0;i=7;i+)ui=hi/(hi+hi+1); /计算三弯矩方程组中的系数ufor(i=1;i=8;i+)vi=1-ui;di=6*(fi+1-fi)/(Xi+1-Xi)-(fi-fi-1)/(Xi-Xi-1)/(hi-1+hi); /计

8、算三弯矩方程组中右边的常数项d r0=2; /从此处开始用追赶法解三对角方程组for(i=1;i=9;i+)li=ui/ri-1;ri=2-li*vi-1;z0=d0; for(i=1;i=0;i-)mi=(zi-vi*mi+1)/2;/此处为前面用追赶法解方程组求得的结果 i=(int)x-1;S=pow(Xi+1-x,3)/(6*hi+1)*mi+pow(x-Xi,3)/(6*hi+1)*mi+1+(fi-(mi*pow(hi+1,2)/6)*(Xi+1-x)/hi+1+(fi+1-(mi+1*pow(hi+1,2)/6)*(x-Xi)/hi+1;/计算三次样条插值函数s（x）s=-pow

9、(Xi+1-x,2)*mi/(2*hi+1)+pow(x-Xi,2)*mi+1/(2*hi+1)-(fi-mi*pow(hi+1,2)/6)/hi+1+(fi+1-mi+1*pow(hi+1,2)/6)/hi+1;/计算s（x）的导数printf(f(4.536)=%fn,S); /输出f(4.536)的结果 printf(f(4.536)=%fn,s); /输出f(4.536)的结果实验结果（截屏图）（3）Matlab语言function Q=san(ssss,p)Q=zeros(2,1);x=1;2;3;4;5;6;7;8;9;10;y=0;0.69314718;1.0986123;1.3

10、862944;1.6094378;1.7917595;1.9459101;2.079445;2.1972246;2.3025851;h=zeros(10,1);d=zeros(10,1);u=zeros(10,1);v=zeros(10,1);r=zeros(10,1);l=zeros(10,1);z=zeros(10,1);m=zeros(10,1);for t=1:1:9;h(t)=x(t+1)-x(t);endd(1)=6/h(1)*(y(2)-y(1)/h(1)-1);d(10)=6/h(9)*(0.1-(y(10)-y(9)/h(9);for t=1:1:8u(t+1)=h(t)/(

11、h(t)+h(t+1);v(t+1)=1-u(t+1);d(t+1)=6/(h(t)+h(t+1)*(y(t+2)-y(t+1)/(x(t+2)-x(t+1)-(y(t+1)-y(t)/(x(t+1)-x(t);endu(10)=1;v(1)=1;r(1)=d(1);for t=2:1:10l(t)=u(t)/r(t-1);r(t)=d(t)-l(t)*v(t-1);endz(1)=d(1);for t=2:1:10z(t)=d(t)-l(t)*z(t-1);endm(10)=z(10)/r(10);for t=9:-1:1m(t)=(z(t)-v(t)*m(t+1)/r(t);endfor

12、t=1:1:10if p=t&p(t+1)Q(:,1)=feval(ssss,p,t,x,m,h,y);breakendendfunction Q=ssss(p,t,x,m,h,y)Q=zeros(2,1);Q(1,1)=(power(x(t+1)-p),3)*m(t)+power(p-x(t),3)*m(t+1)/6+(y(t)-m(t)*h(t)*h(t)/6)*(x(t+1)-p)+(y(t+1)-m(t+1)*h(t)*h(t)/6)*(p-x(t)/h(t);Q(2,1)=(-(power(x(t+1)-p),2)*m(t)+power(p-x(t),2)*m(t+1)/2+(y(t

13、)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6)/h(t);End实验结果（截屏图）（4）问题讨论与分析：样条插值效果比Lagrange插值好,近似误差较小，光滑性较好。本题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。其实编程时对于数组的大小可以弄大点，没必要仔细考虑各数组变量的精确元素个数。同时，编程时也不需要求程序中计算式与原基本公式下标完全一致（例如本题中用追赶法解三对角方程组时计算式中的变量名和书上就不一致），只须数组元素的对应位置（即带入地址）一致即可。第三题 用Romberg算法求（允许误差=0.000 01）。解：（1）

14、解题方法及理论依据：所谓Romberg算法就是采用公式、 按照如下图所示的加工流程计算高精度的romberg积分值。 逐行计算，算完前五行后得，若(准许误差），就把作为积分近似值，否则再算第六行得，若,就把作为积分近似值，否则.,直到时停止计算，并把作为积分近似值。（2） 计算机程序C语言#include #include void main()float f(float x); /函数声明float a=1,b=3,h,m;int n,k;float R10,S10,T10,C10;for(n=1;n=10;n+)m=0; h=(b-a)/(pow(2,(n-1); /求步长h for(k=

15、0;k(pow(2,(n-1);k+) m=m+f(a+(k+0.5)*h); T0=(b-a)*(f(a)+f(b)/2; Tn=0.5*(Tn-1+h*m);/求T Sn-1=(4*Tn-Tn-1)/3;/求s Cn-2=(16*Sn-1-Sn-2)/15; /求C Rn-3=(64*Cn-2-Cn-3)/63;/求Rif(fabs(Rn-3-Rn-4)6 if abs(t(k,4)-t(k-1,4)=15 disp(不收敛); end 调用过程如下： f=inline(3.x).*(x.1.4).*(5.*x+7).*sin(x.2),x) y=rbg(f,1,3,10(-5)试验结果（

16、截屏图）（4）问题讨论与分析：Romberg算法作为一种求积方法加速的这么一种方法，在复化梯形公式的基础上，精确度提高了并且收敛速度加快了。本题中求公式中的 这一部分，其实也可以像被积函数那样做成一个函数来实现。Romberg算法的缺点是：对函数的光滑性要求较高，在计算新分点的值时,这些数值的个数成倍增加。第四题 用定步长四阶Runge-Kutta法求解：h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3)解：（1）解题方法及理论依据： 先不按Taylor公式展开，而是先写成tn处附近的值的线性组合（有待定常数）再按Taylor公

17、式展开，然后确定待定常数，这就是Runge-Kutta法的思想方法。定步长四阶Runge-Kutta公式的古典形式：（2） 计算机程序（C语言）：#include #include float f1(float x,float y,float z) /函数定义float f;f=1;return (f);float f2(float x,float y,float z)/函数定义float f;f=z;return (f);float f3(float x,float y,float z)/函数定义float f;f=1000-1000*y-100*z;return (f);void main

18、()float h=0.0005,k1,k2,k3,k4,n1,n2,n3,n4,m1,m2,m3,m4,x1=0,x2,y1=0,y2,z1=0,z2,t=0; int i=0; /Runge-Kutta算法do t=h*i;if(i=1)x1=x2,y1=y2,z1=z2;/四阶古典Runge-Kutta公式k1=h*f1(x1,y1,z1), m1=h*f2(x1,y1,z1),n1=h*f3(x1,y1,z1),k2=h*f1(x1+k1/2,y1+m1/2,z1+n1/2),m2=h*f2(x1+k1/2,y1+m1/2,z1+n1/2),n2=h*f3(x1+k1/2,y1+m1/

19、2,z1+n1/2),k3=h*f1(x1+k2/2,y1+m2/2,z1+n2/2),m3=h*f2(x1+k2/2,y1+m2/2,z1+n2/2),n3=h*f3(x1+k2/2,y1+m2/2,z1+n2/2),k4=h*f1(x1+k3,y1+m3,z1+n3),m4=h*f2(x1+k3,y1+m3,z1+n3),n4=h*f3(x1+k3,y1+m3,z1+n3),x2=x1+(k1+2*k2+2*k3+k4)/6,y2=y1+(m1+2*m2+2*m3+m4)/6,z2=z1+(n1+2*n2+2*n3+n4)/6;/输出结果if(i=49)printf(y1(%g)=%f y

20、2(%g)=%f y3(%g)=%fn,h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);if(i=89)printf(y1(%g)=%f y2(%g)=%f y3(%g)=%fn,h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);if(i=169)printf(y1(%g)=%f y2(%g)=%f y3(%g)=%fn,h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);if(i=199)printf(y1(%g)=%f y2(%g)=%f y3(%g)=%fn,h*(i+1),x2,h*(i+1),y2,h*(i+1),z2);i+;w

21、hile(i=200);实验结果（截屏图）（3）Maltab语言unction x=mgauss2(A,b,flag)if nargink A(k p,:)=A(p k,:); b(k p,:)=b(p k,:); end m=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag=0,Ab=A,b,endendx=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1 x(k)

22、=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);end format longA=12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743;2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124;-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.0

23、10103;1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585;-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317;0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417;1.742382,-0.784165,-1.056781,2.101023,-3.

24、111223,-0.103458,14.713846,3.123789,-2.213474;3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782;-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001;b=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.

25、6784392;x=mgauss2(A,b);x=x实验结果（截屏图）(4) 问题讨论与分析：Runge-Kutta方法的优点：精度高，不必用别的方法求开始几点的函数值。程序简单，存储量少。方法稳定第五题用列主元素消去法求解AX=b解 (1)解题方法及理论依据：列主元素消去法是建立在Gauss消去法基础上的，Gauss消去法是在较强条件下进行的，即要求Ai非奇异(i=1,2,n-1)。Gauss消去法是按自然顺序消去法，这样的自然顺序选取的主元的解法在计算机中往往误差过大。经过长期经验的积累，我们选取绝对值尽量大的与元素作为主元素，对列选主元的方法称为列主元素消去法。具体做法为：(1)选列主元

26、；(2)交换第一行与列主元所在行；(3)消元计算；(4)回代求解。(2)计算机程序C语言#include#include#define N 9void main()double Fx(double aNN+1,double xN); /*变量及调用函数定义*/int i;double d,xN;double aNN+1=12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743,2.1874369, 2.115237,19.141823,-3.125432,-1.012345,2.189

27、736,1.563849,-0.784165,1.112348,3.123124,33.992318,-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103,-25.173417,1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585,0.84671695,-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,

28、-3.111223,2.121314,1.784317,1.784317,0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417,-86.612343,1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.7138465,3.123789,-2.213474,1.1101230,3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719

29、334,4.446782,4.719345,-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001,-5.6784392; d=Fx(a,x); if(d!=0)for(i=0;iN;i+)printf(x%d=%fn,i+1,xi);double Fx(double aNN+1,double xN) /*用列主元素消去法求解线性方程组的解*/ int i,j,k,r;double b,l,tN;for (k=0;kN;k+)r=k;for(i=k;ifabs(ark) r=i;

30、/得到主元在第r行if(r!=k) /*对调r，k两行*/for(j=k;j=N;j+)b=akj;akj=arj;arj=b;for(i=k+1;iN;i+) /*计算消元因子*/l=aik/akk;for(j=k;j=0;i-)ti=aiN;for(j=i+1;j0,disp(请注意：因为RA=RB，所以此方程组无解.)returnendif RA=RB if RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.) X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)=

31、 B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:n m= B(k,p)/ B(p,p); B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp(请注意：因为RA=RB RA,RB,n,X = Gauss1()实验结果（截屏图）5.3、问题讨论（1）为了使程序简便化,提高运算速率,我们提前把A和b的数值归结到一

32、个数组中，这方便的程序的快速调取和使用；（2）选取主元素和交换行是程序运行的关键,因此在此步骤应着重考虑。（3）简化程序运算过程，使其在获得想要结果的同时普通化，可以通过此程序实现其他方程组的运算。第六题桥梁工程行车道板的计算：背景：设计一座桥首先要重视总体方案、桥形及布置的合理性。上部结构的构造形式、跨径等被确定后，就要进行桥梁各部构件的详细计算。每一座桥梁设计均涉及承载能力计算，并且满足规范要求，承载能力两种状态包括承载能力极限状态设计、正常使用极限状态设计。承载能力极限状态是指桥涵及其构件达到最大承载能力或出现不适于继续承载的变形或变位的状态，正常使用极限状态是指指桥涵及其构件达到正常使用或耐久性的某项限值得状态。计算过程中所用有限元分析软件，比如Midas Civil、Midas FEA等软件，那这些软件的应用与数值分析、编程等又有极为密切的联系，因此掌握数值分析只是第一步，要学会从数值的角度解决问题。意义：钢筋混凝土桥梁的极限状态是桥梁设计不可忽视的一部分，直接承受车辆轮压的钢筋混凝土桥梁，在构造上与桥墩连接在一起，将活载传