数值计算方法上机答案

1、本科实验报告课程名称：实验项目：实验地点：机房专业班级：采矿 1206学号：2012002896学生姓名：康明月指导教师：2014 年7 月 3日一、求非线性方程的根。1、求 方 程 f (x) = x - cos x = 0在 x0 = 1.5附 近 的 是 根 ， 要 求 精 度 满 足x- xk 10-3 .(牛顿切线法)k +1f=inline(x-cos(x); % f(x)df=inline(1+sin(x); %f(x)n=1;x0=input(x0=);del=input(del=);N=input(N=);fprintf(nkx(k);fprintf(n %2d%f ,0,x

2、0);F0=f(x0); dF0=df(x0);while nNif dF0=0fprintf(导数为0,迭代无法继续进行.);return;endx1=x0-F0/dF0;F1=f(x1);dF1=df(x1);if (abs(x1-x0)del) |abs(F1)del)fprintf(n n 结果：%fn,x1);return;endfprintf(n %2d%f ,n,x1);n=n+1;x0=x1;F0=F1;dF0=dF1;endfprintf(nn % d 次迭代后未达到精度要求.n,N);NewtonIterationx0=1.5del=1e-4N=100kx(k)0 1.50

3、00001 0.7844722 0.739519结果：0.7390852、求方程 f (x) = x3 - x2 - 0.8 = 0 在 x0 = 1 附近的是根，求出具有思维有效数字的根近似值.(简单迭代法)clearclcphi=inline(0.8+x2)(1/3); %迭代函数x0=input(x0=);del=input(del=);N=input(N=); n=1;fprintf(n %2d%f,0,x0);while nNx=phi(x0);if abs(x-x0)1s(m)=a(m,n);j=p;while (j2) & (j=m+1) & (jmax1max1=abs(A(i

4、,k);r=i;endendif max1kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);if abs(A(n,n) Gauss(A,b) A = 17.031, -0.615, -2.991,1.007, -1.006,0-1.0, 34.211,-

5、1.0,-2.1,6.3,-1.70,0.5,13.0,-0.5,1.0,-1.54.501,3.11, -3.907, -61.705,12.17, 8.9990.101, -0.812, -0.017,-0.91,4.918,0.11.0,2.0,3.0,4.5,5.0,21.8ans =1.0000-2.00002.9999-4.00015.0000-6.000817 143-1 23- 771 2 10 -17- 2 1143- 4-1 1 - 82- 52-11-11241-11134-1b =- 373、 A = 1317-15 1- 24- 617-1212-1852- 2 13

6、45128-19 2- 735111-1 1- 7 1021精确解x = (4.1992, 0.4955 - 2.0308 7.7589 8.4767 2.7424 7.0443 4.8823)T 其中E = 1E - 4clearclcn=input(n=);A=input(A=);b=input(b=);x=input(x=);epsilon=input(n 精度=);N=input(n 最大迭代次数N=);fprintf(n %d:,0);for i=1:nfprintf(%f,x(i);end%以下是迭代过程for k=1:N%这是第k步迭代，迭代前的向量在x0中，迭代后的向量在x中；

7、 normal=0;for i=1:nt=x(i);x(i)=b(i);for j=1:nif j=ix(i)=x(i)-A(i,j)*x(j);endendx(i)=x(i)/A(i,i);temp=abs(x(i)-t);% 求范数于迭代在同一个循环中； if tempnormalnormal=temp; %这里用的是无穷范数endend %第i不迭代结束；fprintf(n %d: ,k);for i=1:nfprintf(%f,x(i);%输出迭代过程endif normal xx =4.19920.4955-2.03087.75898.47672.74247.04444.882411

8、112341、A =3610411410205153511 5 15 ( pascal矩阵)35701 0b = 0 00精确解x = (5,-10,10,-5,1)T , E = 1E - 6.function x,det,flag = Gauss(A,b)n,m=size(A);nb=length(b);if n=merror(the rows and cloums of A must be equal!);return;endif m=nberror(the rows and cloums of A must be equal the length of b!);return;endfl

9、ag=OK;det=1;x=zeros(n,1);for k=1:(n-1)max1=0;for i=k:nif abs(A(i,k)max1max1=abs(A(i,k);r=i;endendif max1kfor j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:nm=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n

10、);if abs(A(n,n) Gauss(A,b) ans =5-1010-51三、用追赶法求解下列方程组（6 位有效数字）。 4-1000x100 010-1 4 -1 0x2、0-1 4 -1 0x3=0100 -1 40-1x4000-14x5200- 21 x- 0.51- 21 1-1.5 x2 、 =2 1- 21 x9 -1.51- 2 x- 0.5 10x = (6.5,12.5,17,20,21.5,21.5,20,17,12.5,6.5)T解1clearclcn=input(n=);a=input(a=);b=input(b=);c=input(c=);% a(1)=0;

11、c(1)=0; x=input(x=); s=zeros(n,1); t=s;temp=0; for k=1:ns(k)=b(k)-a(k)*temp;% temp t(k)=c(k)/s(k); temp=t(k);end temp=0; for k=1:nx(k)=(x(k)-a(k)*temp)/s(k); temp=x(k);endfor k=n-1:-1:1 x(k)=x(k)-t(k)*x(k+1);endfprintf(nx=n); for k=1:nfprintf(%fn,x(k);27.05138.20513精确解x = 5.7692314.8718endn=5a=0 -1

12、-1 -1 -1b=4 4 4 4 4c=-1 -1 -1 -1 0x=100 0 0 0 200x=27.0512828.2051285.76923114.87179553.717949解2clearclcn=input(n=);a=input(a=);b=input(b=);c=input(c=);% a(1)=0;c(1)=0; x=input(x=); s=zeros(n,1); t=s;temp=0; for k=1:ns(k)=b(k)-a(k)*temp;% temp t(k)=c(k)/s(k); temp=t(k);end temp=0; for k=1:nx(k)=(x(k

13、)-a(k)*temp)/s(k); temp=x(k);endfor k=n-1:-1:1 x(k)=x(k)-t(k)*x(k+1);endfprintf(nx=n); for k=1:nfprintf(%fn,x(k);endn=10a=0 1 1 1 1 1 1 1 1 1 b=-2 -2 -2 -2 -2 -2 -2 -2 -2 -2 c=1 1 1 1 1 1 1 1 1 0x=-0.5 -1.5 -1.5 -1.5 -1.5 -1.5 -1.5 -1.5 -1.5 -0.5x=6.50000012.50000017.00000020.00000021.50000021.50000

14、020.00000017.00000012.5000006.500000四、求解下列积分（E=1E-6）1、12sin xdx 0.6593302、 01 e -x2 dx 0.746824x解 1%复化的梯形公式和复化色Simpson公式计算定积分 f=inline(sin(x)/x); % 函数表达式可以更换； a=input(a=); %积分区间左端点 b=input(b=); %积分区间右端点 n=input(n=); %区间n等分；h=(b-a)/n;temp=f(a);xk=a;for i=1:n-1xk=xk+h;temp=temp+2*f(xk);endtemp=temp+f(

15、b);temp=temp*h/2;fprintf(n复化梯形公式计算结果： %f,temp);fuhuatxa=1b=2n=1000复化梯形公式计算结果： 0.659330解 2%复化的梯形公式和复化色Simpson公式计算定积分 f=inline(exp(-x2); % 函数表达式可以更换； a=input(a=); %积分区间左端点 b=input(b=); %积分区间右端点 n=input(n=); %区间n等分；h=(b-a)/n;temp=f(a);xk=a;for i=1:n-1xk=xk+h;temp=temp+2*f(xk);endtemp=temp+f(b);temp=tem

16、p*h/2;fprintf(n复化梯形公式计算结果： %f,temp);fuhuatxa=0b=1n=1000复化梯形公式计算结果： 0.746824五、求下列矩阵所有最大和最小特征值和相应特征向量。解 1(幂法)function m,u,flag = Pow( A,delta,max1)if nargin3max1=100endif nargin2delta=1e-5;endn=length(A);u=ones(n,1);flag=fail;k=0;m1=0; while k=max1;v=A*u;vmax,i=max(abs(v);m=v(i);u=v/m;ifabs(m-m1) m,u,

17、flag=Pow(A) max1 =100m =3.1738u =0.13090.30710.58990.89801.00000.6815flag =OK(反幂法)function m,u,flag = Pow_in( A,delta,max1)if nargin3max1=100;endif nargin2delta=1e-5;endn=length(A);u=ones(n,1);flag=fail;k=0;m1=0; invA=inv(A);while k=max1v=invA*u;vmax,i=max(abs(v);m=v(i);u=v/m;if abs(m-m1)deltaflag=O

18、K;break;endm1=m;k=k+1;endm=1/m;m,u,flag=Pow_in(A)m =0.8263u =0.1307-0.30680.5896-0.89781.0000-0.6816flag =OK解 2(幂法)function m,u,flag = Pow( A,delta,max1)if nargin3max1=100endif nargin2delta=1e-5;endn=length(A);u=ones(n,1);flag=fail;k=0;m1=0; while k=max1;v=A*u;vmax,i=max(abs(v);m=v(i);u=v/m;ifabs(m-

19、m1) m,u,flag=Pow(A)max1 =100m =-3.6825u =0.5463-0.91911.0000-0.76340.28460.2846-0.76341.0000-0.91910.5463flag =OK(反幂法)function m,u,flag = Pow_in( A,delta,max1)if nargin3max1=100;endif nargin2delta=1e-5;endn=length(A);u=ones(n,1);flag=fail;k=0;m1=0; invA=inv(A);while k=max1v=invA*u;vmax,i=max(abs(v);

20、 m=v(i);u=v/m;if abs(m-m1) yyyy =3.2864Yy 为 3 次插值后的值 y = x + y在区间0,1上八 用欧拉法、预估-校正法和四阶龙格-库塔方法求解初值问题： = 0 y(0)的数值解， h = 0.1 ，并与精确解 y = -x -1+ 2e x 相比较（取 5 位有效数字）。f=inline(y+x,x,y);y=inline(-x-1+2*exp(x);a=input(a=); %区间左端点b=input(b=); %区间右端点h=input(h=); %区间步长；y0=input(y0=);%初始条件；fprintf(n 精确解Euler方法);

21、fprintf( 预估校正方法4阶Runge-Kutta方法);fprintf(n xyyek |yek-y|ymk);fprintf( |ymm-y|yrk |yr(k-y|n);fprintf(%3.1f%8.6f %8.6f %8.6f, a,y0,y0,0);fprintf(%8.6f%8.6f %8.6f %8.6fn,y0,0,y0,0.0);x=a;ye=y0; %Euler 法的初值；ym=y0; % 预估校正方法的初值yr=y0; % 4阶Runge-Kutta方法的初值while (x(b-0.1)ye=ye+h*f(x,ye); %Euler 法在下一个分点的值yp=ym

22、+h*f(x,ym); % 预估校正方法在下一个分点的预报值 ym=ym+(h/2)*(f(x,ym)+f(x+h,yp); % 预估校正方法在下一个分点的值k1=f(x,yr);k2=f(x+h/2,yr+(h/2)*k1);k3=f(x+h/2,yr+(h/2)*k2);k4=f(x+h,yr+h*k3);yr=yr+(h/6)*(k1+2*k2+2*k3+k4);% 4阶Runge-Kutta方法在下一个分点值x=x+h; %下一个分点；yx=y(x);%下一个分点的精确值；fprintf(%3.1f%8.6f %8.6f %8.6f,x,yx,ye,abs(ye-yx);fprintf(%8.6f %8.6f %8.6f %8.6fn,ym,abs(ym-yx),yr,abs(yr-yx); endODEa=0b=1h=0.1y0=1精确