贵州大学-数值分析上机

数值分析上机实验报告学院：土木工程学院专业名称：土木工程姓名：XXXX学号：201XXXXXXX指导老师：XXX老师2017年01月02日第一题一，用Newton法求方程：x7-28x4+14=0在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001.)解：（1）解题方法及理论依据：Newton迭代法定理：设函数在有限区间,ab上二阶导数存在，且满足条件：①()()0fafb；②)(''xf在区间,ab上不变号；③0)('xf；④)()('cfabcf,其中c是,ab中使)(,)(min''bfaf达到的一个。则对任意初始近似值0,xab，由Newton迭代过程1'()(),()kkkkkfxxxxfx0,1,2k所生的迭代序列kx平方收敛于方程()0fx在区间,ab上的唯一解。(2)计算机程序C语言#includestdio.h#includemath.hvoidmain(){intflag=1,c;//设置flag控制while循环（即控制整个程序是否继续运行）while(flag==1){floatx,y,f,f1;inti=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++;//累计迭代的次数printf(x=%f,y=%f\n,x,y);//输出每次迭代的值}printf(\n方程的近似根:y=%f\n,y);printf(\n迭代次数:%d\n,i);printf(是否继续(Y/N)?);//提示是否输入下一个x的初值getchar();c=getchar();if(c=='N'||c=='n')//如果输入N或n，则flag=0退出程序；则继续运行程序flag=0;elseflag=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/fxx;whileabs(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迭代法是一个二阶收敛迭代式，他的几何意义1kx是kx的切线与x轴的交点,故也称为切线法。它是平方收敛的，但它是局部收敛的，即要求初始值与方程的根充分接近，所以在计算过程中需要先确定初始值。本题在讨论中，讨论了区间(0.1,1.9)两端点是否能作为Newton迭代的初值，结果发现0.1不满足条件，而1.9满足，能作为初值0x。第二题二，已知函数值如下表：x12345f(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）解题方法及理论依据：121111121113131)6()6(6)(6)()(iiiiiiiiIiiiiiiihxxhMfhxxhMfMhxxMhxxxS公式112111121112121'1)6(1)6(2)(2)()(iiiiiiiiiiiiiihhMfhhMfMhxxMhxxxS公式2三弯矩方程：nnniiiiiidMMdMMMdMM222111010其矩阵形式：nndddMMM10101-n1n11212212其中),(6),(6),(),(611010101111nnnnniiiiiiiiiiiiixxffhdfxxfhdxxfxxfhhdhhh对于上述矩阵方程，由于其系数矩阵为三对角矩阵，所以用追赶法来解。解三对角方程组的追赶法：设线性方程组Ax=b的系数矩阵A为三对角矩阵nnndacac00dA1211,这时A可以作如下LU分解：nnnrcrcrllLUA001010112112(1)LU分解：首先11dr，对i=2，3，...,n,计算1iiiral,1iiiicldr(2)解Lz=b：首先11zb,对i=2，3，...,n,计算1iiiiylbz(3)解Ux=z：首先nnnrzx,对i=n-1,n-2，....,2,1,计算iiiiirxczx1最后把追赶法求得的唯一解nMMM......,,10，代入上面的公式1和公式2，即可求得xSxS'和。（2）计算机程序C语言#includestdio.h#includemath.hvoidmain(){floatX[10]={1,2,3,4,5,6,7,8,9,10};//输入节点floatf[10]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.94594101,2.079445,2.1972246,2.3025851};//输入节点处的函数值floath[10];floatu[9],v[9],c[9];floatr[10],z[10],l[10],a[10],d[10],m[10];inti;floatx=4.563;floatS,s;for(i=0;i=8;i++)h[i]=X[i+1]-X[i];//计算步长hd[0]=6*((f[0]-f[1])/(X[0]-X[1])-1)/h[1];d[8]=6*(0.1-(f[8]-f[9])/(X[8]-X[9]))/h[8];u[8]=1;v[0]=1;for(i=0;i=7;i++)u[i]=h[i]/(h[i]+h[i+1]);//计算三弯矩方程组中的系数ufor(i=1;i=8;i++){v[i]=1-u[i];d[i]=6*((f[i+1]-f[i])/(X[i+1]-X[i])-(f[i]-f[i-1])/(X[i]-X[i-1]))/(h[i-1]+h[i]);//计算三弯矩方程组中右边的常数项d}r[0]=2;//从此处开始用追赶法解三对角方程组for(i=1;i=9;i++){l[i]=u[i]/r[i-1];r[i]=2-l[i]*v[i-1];}z[0]=d[0];for(i=1;i=9;i++)z[i]=d[i]-l[i]*z[i-1];m[9]=z[9]/2;for(i=8;i=0;i--)m[i]=(z[i]-v[i]*m[i+1])/2;//此处为前面用追赶法解方程组求得的结果i=(int)x-1;S=pow(X[i+1]-x,3)/(6*h[i+1])*m[i]+pow(x-X[i],3)/(6*h[i+1])*m[i+1]+(f[i]-(m[i]*pow(h[i+1],2))/6)*(X[i+1]-x)/h[i+1]+(f[i+1]-(m[i+1]*pow(h[i+1],2)/6))*(x-X[i])/h[i+1];//计算三次样条插值函数s（x）s=-pow(X[i+1]-x,2)*m[i]/(2*h[i+1])+pow(x-X[i],2)*m[i+1]/(2*h[i+1])-(f[i]-m[i]*pow(h[i+1],2)/6)/h[i+1]+(f[i+1]-m[i+1]*pow(h[i+1],2)/6)/h[i+1];//计算s（x）的导数printf(f(4.536)=%f\n,S);//输出f(4.536)的结果printf(f'(4.536)=%f\n,s);//输出f'(4.536)的结果}实验结果（截屏图）（3）Matlab语言functionQ=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);fort=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));fort=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);fort=2:1:10l(t)=u(t)/r(t-1);r(t)=d(t)-l(t)*v(t-1);endz(1)=d(1);fort=2:1:10z(t)=d(t)-l(t)*z(t-1);endm(10)=z(10)/r(10);fort=9:-1:1m(t)=(z(t)-v(t)*m(t+1))/r(t);endfort=1:1:10ifp=t&&p(t+1)Q(:,1)=feval(ssss,p,t,x,m,h,y);breakendendfunctionQ=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);End实验结果（截屏图）（4）问题讨论与分析：①样条插值效果比Lagrange插值好,近似误差较小，光滑性较好。②本题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。③其实编程时对于数组的大小可以弄大点，没必要仔细考虑各数组变量的精确元素个数。同时，编程时也不需要求程序中计算式与原基本公式下标完全一致（例如本题中用追赶法解三对角方程组时计算式中的变量名和书上就不一致），只须数组元素的对应位置（即带入地址）一致即可。第三题用Romberg算法求dxxxxx24.131sin753（允许误差ε=0.00001）。解：（1）解题方法及理论依据：所谓Romberg算法就是采用公式①、②、③、④、⑤按照如下图所示的加工流程计算高精度的romber