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

（数值分析上机实验报告）院 系： 矿业学院 专 业： 矿业工程 班 级： 2015 姓 名： 王 学 号： 2015022 指导教师： 代 第一题1.用Newton法求解方程，在（0.1，1.9）的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。1.1理论依据及方法应用条件Newton迭代法:由一般迭代函数，取s=2时，有，可得二阶迭代序列，此种迭代法称为Newton迭代法。定理:设函数在有限区间a,b上二阶导数存在，且满足条件（）（）在区间a,b上不变号;（）0;（）|/b-a|其中c是a,b中使min|,达到的一个;则对任意时近似值x0a,b,由Newton迭代过程有: k=0,1,2所产生的迭代序列x0平方收敛于方程=0区间a,b上的唯一解。 推论:设函数f(x)满足定理中条件,若选初值，使0,则Newton迭代过程 (k=0,1,2) 产生的迭代序列xk单调收敛于=0的唯一解。1.2计算程序# include # include # include # include using namespace std;double *newton (double a,double b,double eps);/牛顿迭代函数double newtonz (double x);/牛顿迭代子函数void main ()double a=0.1,b=1.9,eps=0.00001,*result;/初始数据coutn 牛顿法解方程：x7-28x4+14=0,在(0.1,1.9)中求近似根，初始值为区间端点，n误差为0.00001。nendl;cout学号：2014021966 姓名：徐林nendl;result = newton (a,b,eps);if (a=result0&result0=b)cout近似根为：result0endl;if (a=result1&result1=b)cout近似根为：result1endl;/-coutn结束,按任意键关闭eps)/代入a迭代计算k+;x2 = x1;x1 = newtonz(x1);/调用牛顿迭代子函数x0 = fabs(x1-x2);re0 = x1;x0 = b, k = 0;while (x0eps)/代入b迭代计算k+;x2 = x1;x1 = newtonz(x1);/调用牛顿迭代子函数x0 = fabs(x1-x2);re1 = x1;return re;1.3计算结果打印1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:M if feval(df,x0)=0 d=2;break else x1=x0-feval(f,x0)/feval(df,x0); end e=abs(x1-x0); x0=x1; if e=eps&abs(feval(f,x1) x0=1.9; eps=0.00001; M=100; x=Newton(f,df,x0,eps,M); vpa(x,7)1.5问题讨论1.需注意的是，要使用Newton迭代法须 满足定理中的条件,，以及0。要用误差范围来控制循环的次数，保证循环的次数和质量，编写程序过程中要注意标点符号的使用，正确运用适当的标点符号，Newton迭代法是局部收敛的，在使用时应先确定初始值，否则所得的解可能不在所要求的范围内。(3)因为newton法求方程是平方收敛的，所以较为精确,但是要求出函数的导数，且必须有二阶导数。第二题2.已知函数值如下表:1234500.693147181.09861231.38629441.60943786789101.79175951.94591012.0794452.19722462.3025851=1=0.1试用三次样条插值求及的近似值。2.1理论依据及方法应用条件 三次样条插值函数可定义为：对于a,b上的一个划分ax0x1x2xn-1=2)如果定义在a,b上的函数S(x),满足（1）.在xi,xi+1上为3次多项式；（2）. S(x), S(x)，S(x)在a,b上连续，则称S(x)在 a,b上划分的3次样条函数，如果对于，还满足，则称为的三次样条插值函数。其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中1,x,x2,x3,(x-xj)+3为基函数,而取B样条函数3(（x-xj）/h)为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X-1,XN+1,则任意三次样条函数可用3(（x-xj）/h)线性组合来表示S(x)= cj3(（x-xj)/h)这样对不同插值问题,若能确定cj由解的唯一性就能求得S(x)。由s(xi)=yi,I=1,2,N s(x0)=y0,s(xN)=yN可得 S(xi)= cj3(（xi-xj）/h)=yi S(x0)=1/hcj3(（x0-xj）/h)=y0 S(xN)=1/hcj3(（xN-xj）/h)=yN2.2计算程序#include#include#define N 10 /*宏定义*/main() 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 dN,uN-1,vN-1,aN-1,cN-1,BN,lN-1,pN,XN; for(j=1;j=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;j=8;j+) vj=hj/(hj-1+hj); for(j=0;j=8;j+) aj=uj; for(j=0;j=8;j+) cj=vj; for(B0=b0,j=1;j=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;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); /*解f(x)的值*/ ds=-X3*pow(x4-t),2)/2/h3+X4*pow(t-x3),2)/2/h3-(y3-X3*h3*h3/6)/h3+(y4-X4*h3*h3/6)/h3; /*解f(x)的值*/ printf(s=%fnds=%fn ,s,ds); /*打印结果*/2.3计算结果打印2.4 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.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=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;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)-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.若要用追赶法求解三对角方程组，三对角阵需要满足：(i=1,2,n)均非奇异，保证有唯一的Doolittle分解;0;2.样条插值效果比Lagrange插值好, 三次样条插值的解存在且唯一，近似误差较小.并且没有Runge现象。第三题3.用Romberg算法求（允许误差=0.00001）。3.1理论依据及方法应用条件 数值积分的Romberg算法计算步骤如下: 当时，就停机。3.2计算程序#include#include#define N 9float f(float x) /*定义函数f(x)*/ float y; y=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(y); main() 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;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)=1e-5) break; i+; printf(The answer is:%f,Ti+10);3.3计算结果打印3.4 MATLAB上机程序function T,n=mromb(f,a,b,eps)if nargineps) J=J+1;h=h/2;S=0; for i=1:n x=a+h*(2*i-1); S=S+feval(f,x); end R(J+1,1)=R(J,1)/2+h*S; for k=1:J R(J+1,k+1)=(4k*R(J+1,k)-R(J,k)/(4k-1); end err=abs(R(J+1,J+1)-R(J+1,J); n=2*n;endR;T=R(J+1,J+1);format longf=(x)(3.x)*(x.1.4)*(5*x+7)*sin(x*x);T,n=mromb(f,1,3,1.e-5)3.5问题讨论1、Romberge算法的优点是:a、把积分化为代数运算,而实际上只需求T1(i),以后用递推可得。b、算法简单且收敛速度快,一般4或5次即能达到要求。c、节省存储量，算出的可存入。2、Romberge算法的缺点是:a、对函数的光滑性要求较高。b、计算新分点的值时,这些数值的个数成倍增加。第四题4.用定步长四阶Runge-Kutta法求解,打印,4.1理论依据及方法应用条件 Runge-Kutta法的基本思想:不是按Taylor公式展开，而是先写成处附近的值的线性组合（有待定系数）,再按Taylor公式展开，然后确定待定常数。 四阶古典Runge-Kutta公式:4.2计算程序#include int main()int i;double h=0.0005;double k1,k2,k3,k4;double y1=0.0,y2=0.0,y3=0.0;for(i=1;i0,这样只需要把上式右端第二项取为2sign(a21)a21 s1即可，从而归纳起来算法步骤为：（2）程序主体#include #include void main( )int i,j,m,r,sign;double A099,s,z,p,n,v,h,l,y9,u9,k,q9,X9,x9=0,0,0,0,0,0,0,0,0,B99,g9;double A99=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.010103,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.741865,0.431637,9.789365,-0.103458,-1.103456,0.238417,1.742382,-0.784165,-1.056781,2.101023,-3.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;double b9=2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392;for(i=0;i9;i+)for(j=0;j9;j+)A0ij=Aij;for(r=0;r7;r+)s=0;for(i=r+1;i0) sign=1;else if(Ar+1r0) sign=-1;for(i=0;i9;i+)if(i=r) ui=0;else if(i=r+1) ui=Ar+1r+sign*s;else ui=Air;for(i=0;i9;i+)yi=0;for(j=0;j9;j+)yi=yi+Aij*uj;yi=yi/l;k=0;for(i=0;i9;i+)k=k+ui*yi;k=k/(2*l);for(i=0;i9;i+)qi=yi-k*ui;for(i=0;i9;i+)for(j=0;j9;j+)Aij=Aij-qi*uj-qj*ui; printf(nHousehold变换矩阵B为:n); for(i=0;i9;i+)printf(n);for(j=0;j9;j+)printf(%f,Aij);（3）结果截屏(4) MATLAB上机程序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)=(b(k)-A(k,k+1:n)*x(k+1:n)/A(k,k);endformat 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.010103;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.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.6784392;x=mgauss2(A,b);x=x（5）问题讨论1、矩阵A是严格对角占优阵，而且还是实数对称阵，满足Householder变换的条件，并且得到的B就是与A相似的三对角阵。2、经过Householder变换矩阵有良好的性质，便于计算分析。