数值分析上机实习报告.docx

1、数值分析上机实习报告（数值分析上机实验报告）院 系： 矿业学院专 业： 矿业工程班 级： 2015姓 名： 王学 号： 2015022指导教师： 代第一题1.用Newton法求解方程X 28x4 14 0，在(0.1 , 1.9 )的近似根(初始近似值取为区间端点, 迭代6次或误差小于0.00001 )。1.1理论依据及方法应用条件Newton迭代法：由一般迭代函数s(x) x ( 1)j ) f (x) f (x)j，取s=2时，有2(x) x丄凶 j 1 j! f (x)可得二阶迭代序列xk 1 xk丄型，此种迭代法称为Newton迭代法。f(xQ定理:设函数在有限区间a,b上二阶导数存在

2、，且满足条件(l)f(a)f(b) 0;(n) f(x)在区间a,b上不变号；(m) f(x)工0;(W) | f(c)|/b-a 0,则Newton迭代过程f (Xk) (k=0,1,2)产生的迭代序列Xk单调收敛于f (x) =0的唯一解af(Xk)1.2计算程序#in elude #in elude #in elude #in elude /牛顿迭代函数牛顿迭代子函数初始数据using n amespaee std;double *n ewt on (double a,double b,double eps); double n ewt onz (double x);void mai n

3、 ()double a=0.1,b=1.9,eps=0.00001,*result;eoutn 牛顿法解方程：xA7-28xA4+14=0,在(0.1,1.9)中求近似根，初始值为区间端点， n误差为0.00001。ne ndl;eout学号：2014021966 姓名：徐林 nendl;result = newton (a,b,eps);if (a=result0&result0=b)cout 近似根为： result0endl;if (a=result1&result1=b)cout 近似根为： result1endl;/ coutn 结束 ,按任意键关闭 eps)k+;x2 = x1;x

4、1 = newtonz(x1); x0 = fabs(x1-x2);re0 = x1;/ 代入 a 迭代计算/调用牛顿迭代子函数x0 = b, k = 0;while (x0eps)k+;x2 = x1;x1 = newtonz(x1);x0 = fabs(x1-x2);re1 = x1;/代入 b 迭代计算/调用牛顿迭代子函数return re;1.3计算结果打印1 C:U5er5Admini Strato Desllop、.D 亡 bug胡亚梅 eyeli.845497Piess anj/ ke to continue1.4MATLAB上机程序function y=Newton(f,df,

5、xO,eps,M) d=0;for k=1:Mif feval(df,xO)=Od=2;breakelsex1= xO-feval(f,xO)/feval(df,xO);ende=abs(x1-x0);x0=x1;if e x0=1.9; eps=0.00001; M=100; x=Newto n(f,df,x0,eps,M); vpa(x,7)1.5问题讨论1.需注意的是，要使用 Newton迭代法须f (x) x7 28x4 14满足定理中的条件I , n ,川，以及f(xc) f(x)0。要用误差范围来控制循环的次数，保证循环的次数和质量，编写程序过程中要注意标点符号的使用，正确运用适当

6、的标点符号， Newton迭代法是局部收敛的，在使用时应先确定初始值，否则所得的解可能不在所要求的范围内(3)因为newton法求方程是平方收敛的，所以较为精确，但是要求出函数的导数，且必须有二阶导数。第二题2.已知函数值如下表：X12345f(x)00.69314718 :1.09861231.38629441.6094378X678910f(x)1.79175951.94591012.0794452.19722462.3025851f(x)f(1)=1f(10)=0.1试用三次样条插值求f (4.563)及f(4.563)的近似值2.1理论依据及方法应用条件 三次样条插值函数可定义为：对于

7、a,b上的一个划分naxoxi X2. Xn-i =2)如果定义在a,b上的函数S(x),满足(1).在x i,x i+i 上为3次多项式；(2) . S(x), S (x),S (x)在a,b上连续，则称S(x)在a,b上划分 的3次样条函数，如果对于f (x) a,b,s(x) f (x)还满足s(xj f (x), i 0n，则称s(x)为f (x)的三次样条插值函数。其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,x 2,x3,(x-x j)+3为基函数, 而取B样条函数Q 3( (x-Xj)/h)为基函数.由于三次样条函数空间是N+3隹的,故我们把分点扩大到X1,X

8、n+1 则任意三次样条函数可用Q 3( (X-Xj) /h)线性组合来表示S(x)= N 1 Cj Q 3( (x-x j)/h)这样对不同插值j 1问题,若能确定Cj由解的唯一性就能求得S(x)。由 s(x i)=y i ,I=1,2,N s (x 0)=y 0 ,s (x N)=y n 可得N 1S(X i )= Cj Q 3( ( Xi -X j )/h)=y i(x)=1/hCj Q 3(X0-Xj )/h)=y(x N)=1/hCj Q 3(XN-Xj )/h)=yd0d1dNd0A yo)0 0d0y y .(j 1 j h 4 h .(j 1 jh .j(yNhN 1yNyjhj

9、 1yj 1)6f(Xj 1，Xj,Xj 1)/*宏定义*/*追赶法求解三弯矩方程*/hj 1hj2.2计算程序#i nclude#in clude#defi ne N 10mai n()float s,ds,t;float dy0=1,dy9=0.1;int j;int xN=1,2,3,4,5,6,7,8,9,10;float yN=0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1972246,2.3025851;int bN=2,2,2,2,2,2,2,2,2,2,hN-1;float

10、dN,uN-1,vN-1,aN-1,cN-1,BN,lN-1,pN,XN;for(j=1;jv=9;j+)hj-1=xj-xj-1;d0=6/h0*(y1/h0-y0/h0-dy0);d9=6/h8*(dy9-y9/h8+y8/h8);for(j=1;j=8;j+) dj=6/(hj-1+hj)*(yj+1/hj-yj/hj-yj/hj-1+yj-1/hj-1);for(u8=1,j=0;j=7;j+)uj=hj-1/(hj-1+hj);for(v0=1,j=1;jv=8;j+)vj=hj/(hj-1+hj);for(j=0;j=8;j+)aj=uj;for(j=0;j=8;j+)cj=vj;

11、for(B0=b0,j=1;jv=9;j+)Bj=bj-aj/Bj-1*cj-1;for(j=1;j=9;j+)lj=aj/Bj-1;for(j=1;j=0;j-)Xj=pj/Bj-cj*Xj+1/Bj;t=4.563;/* 解 f(x) 的值/* 解 f (x)/* 打印结果 */s=X3*pow(x4-t),3)/6/h3+X4*pow(t-x3),3)/6/h3+ (y3-X3*h3*h3/6)*(x4/h3-t/h3)+ (y4-X4*h3*h3/6)*(t/h3-x3/h3);*/ds=-X3*pow(x4-t),2)/2/h3+X4*pow(t-x3),2)/2/h3- (y3-X

12、3*h3*h3/6)/h3+(y4-X4*h3*h3/6)/h3; 的值 */printf(s=%fnds=%fn ,s,ds); 2.3 计算结果打印2.4MATLAB上机程序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.3862944;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=

13、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)/(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;

14、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 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)

15、-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)-m(t)*h(t)*h(t)/6)+(y(t+1)-m(t+1)*h(t)*h(t)/6)/h(t); end2.5 问题讨论1.若要用追赶法求解三对角方程组， 三对角阵需要满足：A(i=1,2,n)均非奇异，保证A有唯一的Doolittle

16、 分解；aG 丰 0;2.样条插值效果比 Lagrange 插值好 , 三次样条插值的解存在且唯一，近似误差较小 . 并且没有 Runge 现象。第二题3.用 Romberg法求：3xx1.4（5x 7）sinx2dx （允许误差&3.1理论依据及方法应用条件弹b ? f（a）f（b） 数值积分的Romberg算法计算步骤如下： 1 （I 1）T1 T124mTmk） Tm=0.00001 )。(0) (0)11 1时,(k 1)Tm 1就停机b a(k 1)21 1f ai 1(2im 1,2, , l k 1,2, , l m 13.2计算程序/*定义函数f（x）*/#i nclude #

17、in clude #defi ne N 9 float f(float x)float y;y=pow(3,x)*pow(x,1.4)*(5*x+7)*si n(x*x);return(y);mai n()float TN+1N+1,hN+1,a=1,b=3,mN+1;int i,l;T10=(b-a)*(f(a)+f(b)/2;l=1;while(l=N)ml=0;for(i=1;i=(pow(2,l-1);i+)ml+=f(a+(2*i-1)*(b-a)/pow(2,l);T1l=(T1l-1+(b-a)* ml/pow(2,l-1)/2;l+;i=1;while(i=N)for(l=1;

18、l=N-i+1;l+)Ti+1l-1=(pow(4,i)*Til-Til-1)/(pow(4,i)-1);hi=Ti0-Ti+10;if(fabs(hi)eps)J=J+1;h=h/2;S=0;for i=1: nx=a+h*(2*i-1);S=S+feval(f,x);endR(J+1,1)=R(J,1)/2+h*S;for k=1:JR(J+1,k+1)=(4Ak*R(J+1,k)-R(J,k)/(4Ak-1);enderr=abs(R(J+1,J+1)-R(J+1,J);n=2*n;endR;T=R(J+1,J+1);format longf=(x)(3.Ax)*(x.A1.4)*(5*

19、x+7)*si n(x*x);T, n=mromb(f,1,3,1.e-5)3.5问题讨论1、 Romberge算法的优点是：a、 把积分化为代数运算，而实际上只需求,以后用递推可得b、 算法简单且收敛速度快，一般4或5次即能达到要求。c、 节省存储量，算出的可存入。2、 Romberge算法的缺点是：a、 对函数的光滑性要求较高。b、 计算新分点的值时，这些数值的个数成倍增加。第四题4.用定步长四阶Runge-Kutta法求解厂 dy,/dt 1dyz/dt y3dy3/ dt 1000 1000 y2 100 y3yi(0) 0y2(0) 0I y3(0) 0h 0.0005,打印 yi

20、(0.025) , y (0.045) , % (0.085) , %(0.1), (i 1,2,3)4.1理论依据及方法应用条件Runge-Kutta法的基本思想：旳1不是按Taylor公式展开，而是先写成tn处附近的值的线性组合（有 待定系数）,再按Taylor公式展开，然后确定待定常数。四阶古典Runge-Kutta公式：Yn 1Yn (K1 2K62 2K3 K4)K1hF(Xn,Yn)K2hF(Xnh,Yn22K1)2K3hF(Xnh,Yn2丄心）2K4hF(Xnh,Yn21-K3)24.2计算程序#include int mai n()int i;double h=0.0005;d

21、ouble k1,k2,k3,k4;double y1=0.0,y2=0.0,y3=0.0;for(i=1;i=200;i+)k1= k2=k3=k4=h*1.0;y1+=(k1+2*k2+2*k3+k4)/6;k仁 k2=k3=k4=h*y3; y2+=(k1+2*k2+2*k3+k4)/6;k仁h*(1000-1000*y2-100*y3); k2=h*(1000-1000*y2-100*(y3+0.5*k1); k3=h*(1000-1000*y2-100*(y3+0.5*k2); k4=h*(1000-1000*y2-100*(y3+k3); y3+=(k1+2*k2+2*k3+k4)

22、/6;if(i=50)prin tf(ny1(0.025)=%f con ti nue;if(i=90)prin tf(ny1(0.045)=%f con ti nue;if(i=170)prin tf(ny1(0.085)=%f con ti nue;if(i=200)prin tf(ny1(0.100)=%fy2(0.025)=%fy2(0.045)=%fy2(0.085)=%fy2(0.100)=%fy3(0.025)=%f,y1,y2,y3);y3(0.045)=%f,y1,y2,y3);y3(0.085)=%f,y1,y2,y3);y3(0.100)=%fnn,y1,y2,y3);4

23、.3计算结果打印4.4 MATLAB上机程序function Y=R_K(df1,a,b,h)m=(b-a)/h;Y=zeros(3,1);S=zeros(3,1);K=zeros(3,4);x=a;y1=a;y2=a;y3=a;for n=1:mK(:,1)=feval(df1,x,y1,y2,y3);x=x+0.5*h;S(:,1)=Y+0.5*h.*K(:,1);y1=S(1,1);y2=S(2,1);y3=S(3,1);K(:,2)=feval(df1,x,y1,y2,y3);S(:,1)=Y+0.5*h.*K(:,2);y1=S(1,1);y2=S(2,1);y3=S(3,1);K(

24、:,3)=feval(df1,x,y1,y2,y3);x=x+0.5*h;S(:,1)=Y+h.*K(:,3);y1=S(1,1);y2=S(2,1);y3=S(3,1);K(:,4)=feval(df1,x,y1,y2,y3);Y=Y+h.*(K(:,1)+2.*K(:,2)+2.*K(:,3)+K(:,4)/6;endfunction Z=df1(x,y1,y2,y3)Z=zeros(3,1);Z(1)=0*x+0*y1+0*y2+0*y3+1;Z(2)=0*x+0*y1+0*y2+1*y3;Z(3)=0*x+0*y1-1000*y2-100*y3+1000;end4.5问题讨论1.定步长四阶runge-kutta法稳定，精度高，可根据有y f (t, y)变化的情况与需要的精度自动修改步长，误差小且程序简单，存储量少2.但是Runge-Kutta法需要每步都计算函数值f(t,y)四次,在函数较复杂时，工作量就会变得较大可靠性有待核查。第五题5.已知A与b