數(shù)值分析實(shí)驗(yàn)報(bào)告三

1、數(shù)學(xué)與信息工程學(xué)院實(shí) 驗(yàn) 報(bào) 告課程名稱： 計(jì)算方法 實(shí) 驗(yàn) 室： 7404 實(shí)驗(yàn)臺(tái)號(hào)： 班 級(jí)： 姓 名： 實(shí)驗(yàn)日期： 2014 年 5月 21 日 實(shí)驗(yàn)名稱非線性方程求解及方程組直接解法 實(shí)驗(yàn)?zāi)康暮?要 求1.二分法的Matlab實(shí)現(xiàn)；2.牛頓法的Matlab實(shí)現(xiàn)；3.牛頓下山法、割線法、艾特金加速法、重根迭代法、非線性方程 組牛頓法中任選其一。4.列主元消元法的Matlab實(shí)現(xiàn)；5.二選一：（1）平方根法的Matlab實(shí)現(xiàn)；（2）追趕法的Matlab實(shí)現(xiàn)。 實(shí)驗(yàn)內(nèi)容和步驟：實(shí)驗(yàn)內(nèi)容：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、用迭代法求方程附近的一個(gè)根。方程寫成下 列等價(jià)形式，并建立相應(yīng)的迭代公式。，迭代公式4、用列主元消法解線性方程組Ax=b。5、用平方根法解線性

3、方程組Ax=b。（1）（2）實(shí)驗(yàn)步驟：1） 實(shí)驗(yàn)編程2） 運(yùn)行結(jié)果二分法1） 實(shí)驗(yàn)編程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); %程序中調(diào)用的fun.m 為函數(shù) if ya* yb>0, disp('注意：ya*yb>0,請(qǐng)重新調(diào)整區(qū)間端點(diǎn)a和b.'), returnendmax1=-1+ceil(log(b-a)- log(abtol)/ log(2); % ceil(u)是大于u的最小取整數(shù)for k=1:

4、 max1+1 a;ya=fun(a); 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*yx>0 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

5、=-4:0.1:4;y=exp(x) +10*x-20;plot(x,y)gridk,x,wuca,yx=erfen (1,1,10-5)2） 運(yùn)行結(jié)果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.

6、3831ans = 4.0000 0 0.1250 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

7、.0358 -0.0076ans = 9.0000 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

8、0.0001 -0.0021 0.0006 -0.0008ans = 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.0905wuc

9、a =7.6294e-006yx =-2.5908e-005牛頓法1）實(shí)驗(yàn)編程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=

10、fnq(x(i); (i-1) xk yk piancha xdpianchaif (abs(yk)<ftol)&(piancha<tol)|(xdpiancha< tol) k=i-1; xk=x(i);(i-1) xk yk piancha xdpiancha return;endend if i>gxmax disp('請(qǐng)注意：迭代次數(shù)超過給定的最大值gxmax。') k=i-1; xk=x(i);(i-1) xk yk piancha xdpiancha return;end (i-1),xk,yk,piancha,xdpiancha&#

11、39;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. 在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）

12、 運(yùn)行結(jié)果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 = 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.

13、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 = 0.0056埃特金加速收斂方法1） 實(shí)驗(yàn)編程1. 新建M文件diedail.m輸入：function k,piancha,xdpiancha,xk=diedai1(x0,k)% 輸入的量-x0是初始值,k是迭代次數(shù)x(1)=x0; for i=1:k x(i+

14、1)=fun1(x(i);%程序中調(diào)用的fun1.m為函數(shù)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)&(xdpiancha>0.5)&(k>3) disp('注意：此迭代序列發(fā)散，請(qǐng)重新輸入新的迭代公式') return end if (piancha < 0.001)&(xdpiancha< 0.0000005)&

15、amp;(k>3) 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)+eps);5. 建立名為fun1.m的M文件function f=fun1(x) f=x-(fai

16、(x)-x)2/(fai(fai(x)-2*fai(x)+x);6. 在MATLAB工作窗口輸入：x0=1.5;k=3;diedai1(x0,k)2）運(yùn)行結(jié)果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） 實(shí)驗(yàn)編程新建M文件gaos2.m輸入：function RA,RB,n,X=gaos2(A,b)B=A b; n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA;if

17、zhica>0,disp('請(qǐng)注意：因?yàn)镽A=RB，所以此方程組無解.')returnendif RA=RB if RA=ndisp('請(qǐng)注意：因?yàn)镽A=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); %計(jì)算乘數(shù) B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend b=B(1:n,n+1);A=B(1:n

18、,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('請(qǐng)注意：因?yàn)镽A=RB<n，所以此方程組有無窮多解.')endend2. 在MATLAB工作窗口輸入：A=-3 2 6;10 -7 0;5 -1 5;>> b=4;7;6;>> RA,RB,n,X =gaos2 (A,b) 2）運(yùn)行結(jié)果>> RA,RB,n,X =gaos2 (A,b)請(qǐng)注意：因?yàn)镽A=RB=n，所以此方程組有唯一解.Warnin

19、g: Divide by zero.> In gaos2 at 24RA = 3RB = 3n = 3X = NaN -1 1平方根法1） 實(shí)驗(yàn)編程1.新建M文件pfg.m輸入：function f=pfg(A,b) if A=A' disp('請(qǐng)注意：因?yàn)锳不是對(duì)稱矩陣,所以此方程組無法分解') end D p=chol(A); if p=0 disp('請(qǐng)注意：因?yàn)锳不是正定矩陣,所以此方程組無法分解') endn=length(b);s=0; c=0;d=0;e=0;L=zeros(n,n);y=zeros(1,n);x=zeros(1,n)

20、;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,