数值分析实验报告三

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 r

2、oot of with accuracy .(1)、Use the bisection method starting with the interval 0,1.(2)、Use the iteration method ,the initial value .2、The equation has two roots near 0.1.Determine them by means of Newtons method.(with accuracy )3、用迭代法求方程附近的一个根。方程写成下 列等价形式，并建立相应的迭代公式。，迭代公式4、用列主元消法解线性方程组Ax=b。5、用平方根法解线性

3、方程组Ax=b。（1）（2）实验步骤：1） 实验编程2） 运行结果二分法1） 实验编程1. 新建M文件newtonqx.m输入：function k,x,wuca,yx=erfen(a,b,abtol)a(1)=a; b(1)=b; ya=fun(a(1); yb=fun(b(1); %程序中调用的fun.m 为函数 if ya* yb0, disp(注意：ya*yb0,请重新调整区间端点a和b.), returnendmax1=-1+ceil(log(b-a)- log(abtol)/ log(2); % ceil(u)是大于u的最小取整数for k=1: max1+1 a;ya=fun(a

4、); b;yb=fun(b); x=(a+b)/2; yx=fun(x); wuca=abs(b-a)/2; k=k-1; k,a,b,x,wuca,ya,yb,yx if yx=0 a=x; b=x; elseif yb*yx0 b=x;yb=yx; else a=x; ya=yx; end if b-a abtol , return, endendk=max1; x; wuca; yx=fun(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-2

5、0;plot(x,y)gridk,x,wuca,yx=erfen (1,1,10-5)2） 运行结果ans = 0 -1.0000 1.0000 0 1.0000 -11.6321 10.7183 -1.0000ans = 1.0000 0 1.0000 0.5000 0.5000 -1.0000 10.7183 4.6487ans = 2.0000 0 0.5000 0.2500 0.2500 -1.0000 4.6487 1.7840ans = 3.0000 0 0.2500 0.1250 0.1250 -1.0000 1.7840 0.3831ans = 4.0000 0 0.1250

6、0.0625 0.0625 -1.0000 0.3831 -0.3105ans = 5.0000 0.0625 0.1250 0.0938 0.0313 -0.3105 0.3831 0.0358ans = 6.0000 0.0625 0.0938 0.0781 0.0156 -0.3105 0.0358 -0.1375ans = 7.0000 0.0781 0.0938 0.0859 0.0078 -0.1375 0.0358 -0.0509ans = 8.0000 0.0859 0.0938 0.0898 0.0039 -0.0509 0.0358 -0.0076ans = 9.0000

7、0.0898 0.0938 0.0918 0.0020 -0.0076 0.0358 0.0141ans = 10.0000 0.0898 0.0918 0.0908 0.0010 -0.0076 0.0141 0.0033ans = 11.0000 0.0898 0.0908 0.0903 0.0005 -0.0076 0.0033 -0.0021ans = 12.0000 0.0903 0.0908 0.0906 0.0002 -0.0021 0.0033 0.0006ans = 13.0000 0.0903 0.0906 0.0905 0.0001 -0.0021 0.0006 -0.0

8、008ans = 14.0000 0.0905 0.0906 0.0905 0.0001 -0.0008 0.0006 -0.0001ans = 15.0000 0.0905 0.0906 0.0905 0.0000 -0.0001 0.0006 0.0002ans = 16.0000 0.0905 0.0905 0.0905 0.0000 -0.0001 0.0002 0.0001ans = 17.0000 0.0905 0.0905 0.0905 0.0000 -0.0001 0.0001 -0.0000k =17x =0.0905wuca =7.6294e-006yx =-2.5908e

9、-005牛顿法1）实验编程1. 新建M文件newtonqx.m输入：%输入：初始值x0;近似值xk 的精度tol； f(xk)的精度tol。function k,xk,yk,piancha,xdpiancha=newtonqx(x0,tol,ftol,gxmax)x(1)=x0; for i=1: gxmax x(i+1)=x(i)-fnq(x(i)/(dfnq(x(i)+eps); piancha=abs(x(i+1)-x(i); xdpiancha= piancha/( abs(x(i+1)+eps); i=i+1;xk=x(i);yk=fnq(x(i); (i-1) xk yk pian

10、cha xdpianchaif (abs(yk)ftol)&(pianchatol)|(xdpianchagxmax disp(请注意：迭代次数超过给定的最大值gxmax。) k=i-1; xk=x(i);(i-1) xk yk piancha xdpiancha return;end (i-1),xk,yk,piancha,xdpiancha;2. 建立名为fnq.m的M文件function y=fnq(x) y=2*x.4+24*x.3+61*x.2-16*x+1;3. 建立名为dfnq.m的M文件function y=dfnq(x) y=8*x.3+72*x.2+122*x-16;4.

11、在MATLAB工作窗口输入：x=0.5:0.1:4; y=2*x.4+24*x.3+61*x.2-16*x+1; plot(x,y) grid k,xk,yk,piancha,xdpiancha=newtonqx(0.5,0.001, 0.001,100)2） 运行结果ans = 1.0000 0.3223 3.0037 0.1777 0.5515ans = 2.0000 0.2256 0.7752 0.0967 0.4287ans = 3.0000 0.1748 0.1972 0.0508 0.2903ans = 4.0000 0.1488 0.0497 0.0260 0.1751ans =

12、 5.0000 0.1356 0.0125 0.0132 0.0973ans = 6.0000 0.1289 0.0031 0.0066 0.0514ans = 7.0000 0.1256 0.0008 0.0033 0.0262ans = 8.0000 0.1240 0.0002 0.0016 0.0129ans = 9.0000 0.1233 0.0000 0.0007 0.0056ans = 9.0000 0.1233 0.0000 0.0007 0.0056k = 9xk = 0.1233yk = 3.4012e-005piancha = 6.9657e-004xdpiancha =

13、0.0056埃特金加速收敛方法1） 实验编程1. 新建M文件diedail.m输入：function k,piancha,xdpiancha,xk=diedai1(x0,k)% 输入的量-x0是初始值,k是迭代次数x(1)=x0; for i=1:k x(i+1)=fun1(x(i);%程序中调用的fun1.m为函数y=(x) piancha= abs(x(i+1)-x(i); xdpiancha=piancha/( abs(x(i+1)+eps); i=i+1;xk=x(i);(i-1) piancha xdpiancha xkendif (piancha 1)&(xdpiancha0.5)

14、&(k3) disp(注意：此迭代序列发散，请重新输入新的迭代公式) return end if (piancha 0.001)&(xdpiancha3) disp(祝贺您！此迭代序列收敛，且收敛速度较快) return endp=(i-1) piancha xdpiancha xk;2. 建立名为fnq.m的M文件function f=fnq(x) f=x.3+2*x.2+10*x-20;3. 建立名为dfnq.m的M文件function f=dfnq(x); f=3*x2+4*x+10;4. 建立名为fai.m的M文件function f=fai(x);f=x-fnq(x)/(dfnq(x

15、)+eps);5. 建立名为fun1.m的M文件function f=fun1(x) f=x-(fai(x)-x)2/(fai(fai(x)-2*fai(x)+x);6. 在MATLAB工作窗口输入：x0=1.5;k=3;diedai1(x0,k)2）运行结果ans = 1.0000 0.1314 0.0960 1.3686ans = 2.0000 0.0002 0.0001 1.3688ans = 3.0000 0.0000 0.0000 1.3688ans = 3列主元消元法1） 实验编程新建M文件gaos2.m输入：function RA,RB,n,X=gaos2(A,b)B=A b;

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

17、endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A(n,n); for q=n-1:-1:1 X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q); endelse disp(请注意：因为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 = 3n = 3X = NaN -1 1平方根法1） 实验编程

18、1.新建M文件pfg.m输入：function f=pfg(A,b) if A=A disp(请注意：因为A不是对称矩阵,所以此方程组无法分解) end D p=chol(A); if p=0 disp(请注意：因为A不是正定矩阵,所以此方程组无法分解) endn=length(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:n L(t,1)=A(t,1)/L(1,1);endfor i=3:n j=i-1; for k=1:j-1 s=s+L(j,k)2; c=c+L(i,k)*L(j,k); end L(j,j)=sqrt(A(j,j)-s); L(i,j)=(A(i,j)-c)/L(j,j);ends=0for k=1:n-1 s=s+L(n,k)2; end L(n,n)=sqrt(A(n,n)-s);Ly(1)=b(1)/L(1,1);for m=2:n d=0; for k=1:m-1 d=d+L(m,k)*y(k); end y(m)=(b(m)-d)/L(m,m);endyx(n)=y(n)/L(n,n);for i=n-1:-1:1 e=0 for k=i+1:n e=e+L(