数值分析实验报告三综述

1、数学与信息工程学院实验报告课程名称： 计算方法实验室： 7404实验台号：班 级：姓 名：实验日期： 2014 年 5 月 21 日实验名称非线性方程求解及方程组直接解法实验目的和要求1. 二分法的Matlab实现；2. 牛顿法的Matlab实现；3牛顿下山法、割线法、艾特金加速法、重根迭代法、非线性方程 组牛顿法中任选其一。4. 列主兀消兀法的Matlab头现；5. 二选一：（1）平方根法的Matlab实现；（2）追赶法的Matlab实 现。实验内容和步骤：实验内容：1、 Compare the number of computation for finding the root ofex

2、iox 2 0 with accuracy 105.(1) 、Use the bisection method starting with the interval 0,1.(2) 、 Use the iteration method (2 eXk)/10 ,the initial valueX。0.2、The equation 2x 24x61x16x 1 0 has two roots near 0.1.Determ ine them by means of Newt on method.8(with accuracy | 耳 xk 11 10)3、用迭代法求方程x3 x2 1 0在X0

3、1.5附近的一个根。方程写成下列等价形式，并建立相应的迭代公1式。x 1；，迭代公式Xk 1 112xXk4、用列主元消法解线性方程组Ax=lbo3264A 107 0 ,b 751565、用平方根法解线性方程组Ax=bo67593211(1：)A 7 13 8 ,b 10(2)A 254 ,b2586914474实验步骤：1) 实验编程2) 运行结果二分法1)实验编程1. 新建M文件newtonqx.m输入：fun ctio n k,x,wuca,yx=erfe n(a,b,abtol)a(1)=a; b(1)=b;ya=fun(a(1); yb=fun(b(1);%程序中调用的 fun.m

4、 为函数if ya* yb0,disp(注意：ya*yb0,请重新调整区间端点a和b. ), return endmax仁-1+ceil(log(b-a)- log(abtol)/ log(2);% ceil(u)是大于 u的最小取整数for k=1: max1+1a;ya=fu n( a); b;yb=fu n( b); x=(a+b)/2;yx=fun(x); wuca=abs(b-a)/2; k=k-1;k,a,b,x,wuca,ya,yb,yxif yx=0a=x; b=x;elseif yb*yx0 b=x;yb=yx;elsea=x; ya=yx;endif b-a abtol ,

5、 return , endendk=max1; x; wuca; yx=fu n( x);2. 建立名为fun.m的M文件function fun=fun(x)fun=exp(x)+10*x-2;3在MATLAB 工作窗口输入：x=-4:0.1:4; y=exp(x) +10*x-20;Plot(x,y) grid k,x,wuca,yx=erfen ( - 1,1,10A-5)2)运行结果80d-3-2N1734ans0-1.00001.000001.0000-11.632110.7183-1.0000ans1.000001.00000.50000.5000-1.000010.71834.6

6、487ans2.000000.50000.25000.2500-1.00004.64871.7840ans3.000000.25000.12500.1250-1.00001.78400.3831ans4.000000.12500.06250.0625-1.00000.3831-0.3105ans5.00000.06250.12500.09380.0313-0.31050.38310.0358ans=6.00000.06250.09380.07810.0156-0.31050.0358-0.1375ans =7.00000.07810.09380.08590.0078-0.13750.0358-

7、0.0509ans =8.00000.08590.09380.08980.0039-0.05090.0358-0.0076ans =9.00000.08980.09380.09180.0020-0.00760.03580.0141ans =10.00000.08980.09180.09080.0010-0.00760.01410.0033ans =11.00000.08980.09080.09030.0005-0.00760.0033-0.0021ans =12.00000.09030.09080.09060.0002-0.00210.00330.0006ans =13.00000.09030

8、.09060.09050.0001-0.00210.0006-0.0008ans =14.00000.09050.09060.09050.0001-0.00080.0006-0.0001ans =15.00000.09050.09060.09050.0000-0.00010.00060.0002ans =16.00000.09050.09050.09050.0000-0.00010.00020.0001ans =17.00000.09050.09050.09050.0000-0.00010.0001-0.0000k =17x =0.0905wuca =7.6294e-006yx =-2.590

9、8e-005牛顿法1)实验编程1.新建M文件newtonqx.m输入：%俞入：初始值x0;近似值xk的精度tol ; f(xk)的精度tol 。fun cti onk,xk,yk,pia ncha,xdpia ncha=n ewt on qx(x0,tol,ftol,gxmax)x(1)=x0;for i=1: gxmaxx(i+1)=x(i)-fnq(x(i)/(dfnq(x(i)+eps); pia ncha=abs(x(i+1)_x(i); xdpia ncha= pia ncha/( abs(x(i+1)+eps); i=i+1;xk=x(i);yk=fnq(x(i); (i-1) x

10、k yk pia ncha xdpia nchaif (abs(yk)vftol)&(pianchatol)|(xdpianchagxmaxdisp(请注意：迭代次数超过给定的最大值gxmax。)k=i-1; xk=x(i);(i-1) xk yk pia ncha xdpia nchareturn ;end(i-1),xk,yk,pia ncha,xdpia ncha;2. 建立名为fnq.m的M文件fun cti on y=fnq(x)y=2*x.A4+24*x.A3+61*x.A2-16*x+1;3. 建立名为dfnq.m的M文件fun cti on y=dfnq(x)y=8*x.A3+

11、72*x.A2+122*x-16;4. 在MATLAB 工作窗口输入：x=0.5:0.1:4; y=2*x.A4+24*x.A3+61*x.A2-16*x+1; plot(x,y)grid k,xk,yk,pia ncha,xdpia ncha=n ewto nqx(0.5,0.001,0.001,100)2)运行结果0511525335ans1.00000.32233.00370.17770.5515ans2.00000.22560.77520.09670.4287ans3.00000.17480.19720.05080.2903ans4.00000.14880.04970.02600.17

12、51ans5.00000.13560.01250.01320.0973ans6.00000.12890.00310.00660.0514ans7.00000.12560.00080.00330.0262ans8.00000.12400.00020.00160.0129ans9.00000.12330.00000.00070.0056ans9.00000.12330.00000.00070.00560dUJU10005002血20QCinnk = 9xk = 0.1233yk = 3.4012e-005 pia ncha = 6.9657e-004xdpia ncha = 0.0056埃特金加速

13、收敛方法1)实验编程1. 新建M文件diedail.m 输入：fun cti onk,pia ncha,xdpia ncha,xk=diedai1(x0,k)%输入的量-x0是初始值,k是迭代次数x(1)=x0;for i=1:kx(i+1)=fun1(x(i);%程序中调用的 fun 1.m 为函数 y= $ (x)pia ncha= abs(x(i+1)-x(i); xdpia ncha=pia ncha/( abs(x(i+1)+eps);i=i+1;xk=x(i);(i-1) pia ncha xdpia ncha xk endif (piancha 1)&(xdpiancha0.5)

14、&(k3)disp(注意：此迭代序列发散，请重新输入新的迭代公式)returnendif (piancha 0.001)&(xdpiancha3)disp(祝贺您！此迭代序列收敛，且收敛速度较快)returnendp=(i-1) pia ncha xdpia ncha xk;2. 建立名为fnq.m的M文件fun cti on f=fnq(x)f=x.A3+2*x.A2+10*x-20;3. 建立名为dfnq.m的M文件fun cti on f=dfnq(x);f=3*xA2+4*x+10;4. 建立名为fai.m的M文件fun cti on f=fai(x);f=x-fnq(x)/(dfnq

15、(x)+eps);5. 建立名为fun1.m的M文件fun cti on f=fun 1(x)f=x-(fai(x)-x)A2/(fai(fai(x)-2*fai(x)+x);6.在 MATLAB工作窗口输入：x0=1.5;k=3;diedai1(xO,k)2)运行结果ans =1.00000.13140.09601.3686ans =2.00000.0002 0.00011.3688ans =3.00000.0000 0.00001.3688ans =3列主元消元法1)实验编程新建M文件gaos2.m输入：fun ctionRA,RB, n,X=gaos2(A,b)B=A b; n=le n

16、gth(b); RA=ra nk(A);RB=ra nk(B);zhica=RB-RA;if zhica0,disp(请注意:因为RA=RB，所以此方程组无解.)returnendif RA=RBif RA=ndisp(请注意：因为RA=RB=n，所以此方程组有唯一解.)X=zeros (n ,1);t=1;for p= 1:n-1C1,l1=max(abs(B(:,p);%列主元B(t,l1,:)=B(l1,t,:);t=t+1;for k=p+1: nm= B(k,p)/ B(p,p);%计算乘数B(k,p: n+1)= B(k,p: n+1)_m* B(p,p: n+1);endendb

17、=B(1: n,n+1);A=B(1: n,1: n); X( n)=b( n)/A( n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1: n)*X(q+1: n)/A(q,q);endelsedisp(请注意：因为RA=RB b=4;7;6; RA,RB,n,X =gaos2 (A,b)2)运行结果 RA,RB ,n,X =gaos2 (A,b)请注意：因为RA=RB=n所以此方程组有唯一解.Warning: Divide by zero. In gaos2 at 24RA =3RB =3X =NaN-11平方根法1)实验编程1.新建M文件pfg.m输入：f

18、un cti onf=pfg(A,b)if A=Adisp( 请注意：因为A不是对称矩阵，所以此方程组无法分解)endD p=chol(A);if p=0disp( 请注意：因为A不是正定矩阵，所以此方程组无法分解)endn=len gth(b);s=0; c=0;d=0;e=0;L=zeros( n,n);y=zeros(1, n);x=zeros(1, n);L(1,1)=sqrt(A(1,1);for t=2:nL(t,1)=A(t,1)/L(1,1);endfor i=3:nj=i-1;for k=1:j-1 s=s+L(j,kF2; c=c+L(i,k)*L(j,k);endL(j,j)=sqrt(A(j,j)-s);L(i,j)=(A(i,j)-c)/L(j,j);ends=0for k=1:n-1s=s+L( n,kF2;endL(n,n )=sqrt(A( n,n )-s);Ly(1)=b(1)/L(1,1);for m=2:nd=0;for k=1:m-1 d=d+L(m,k)*y(k); end y(m)=(b(m)-d)/L(m,m); end y x(n )=y( n)/L( n,n); for i=n-1:-1:1