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

1、 數(shù)值分析實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)一、誤差分析誤差問題是數(shù)值分析的基礎(chǔ)，又是數(shù)值分析中一個困難的課題。在實(shí)際計(jì)算中，如果選用了不同的算法，由于舍入誤差的影響，將會得到截然不同的結(jié)果。因此，選取算法時注重分析舍入誤差的影響，在實(shí)際計(jì)算中是十分重要的。同時，由于在數(shù)值求解過程中用有限的過程代替無限的過程會產(chǎn)生截?cái)嗾`差，因此算法的好壞會影響到數(shù)值結(jié)果的精度。一、實(shí)驗(yàn)?zāi)康?、 通過上機(jī)編程，復(fù)習(xí)鞏固以前所學(xué)程序設(shè)計(jì)語言及上機(jī)操作指令；2、 通過上機(jī)計(jì)算，了解誤差、絕對誤差、誤差界、相對誤差界的有關(guān)概念；3、 通過上機(jī)計(jì)算，了解舍入誤差所引起的數(shù)值不穩(wěn)定性。二、實(shí)驗(yàn)任務(wù)對，計(jì)算定積分 .算法1：利用遞推公式 , ,

2、取 .算法2：利用遞推公式 .注意到,取 .思考：從計(jì)算結(jié)果看，哪個算法是不穩(wěn)定的，哪個算法是穩(wěn)定的。算法1：t=log(6.0)-log(5.0);n=0;y=zeros(1,21);y(1)=t;for k=2:21y(k)=1/k-5*y(k-1);n=n+1;endy(1:6)y(7:21)運(yùn)行結(jié)果：ans = 0.1823 -0.4116 2.3914 -11.7069 58.7343 -293.5049算法2：y=zeros(21,1);n=1;y1=(1/105+1/126)/20;for k=21:-1:2y(k-1)=1/(5*k)-y(k)/5;n=n+1;end運(yùn)行結(jié)果：

3、y =0.0884 0.0580 0.0431 0.0343 0.0285 0.0243 0.0212 0.0188 0.0169 0.0154 0.0141 0.0130 0.0120 0.0112 0.0105 0.0099 0.0093 0.0089 0.0081 0.0095 0由數(shù)據(jù)對比可知，算法2較為穩(wěn)定。實(shí)驗(yàn)二、插值法插值法是函數(shù)逼近的一種重要方法，它是數(shù)值積分、微分方程數(shù)值解等數(shù)值計(jì)算的基礎(chǔ)與工具，其中多項(xiàng)式插值是最常用和最基本的方法。拉格朗日插值多項(xiàng)式的優(yōu)點(diǎn)是表達(dá)式簡單明確，形式對稱，便于記憶，它的缺點(diǎn)是如果想要增加插值節(jié)點(diǎn)，公式必須整個改變，這就增加了計(jì)算工作量。而牛頓插值

4、多項(xiàng)式對此做了改進(jìn)，當(dāng)增加一個節(jié)點(diǎn)時只需在原牛頓插值多項(xiàng)式基礎(chǔ)上增加一項(xiàng)，此時原有的項(xiàng)無需改變，從而達(dá)到節(jié)省計(jì)算次數(shù)、節(jié)約存儲單元、應(yīng)用較少節(jié)點(diǎn)達(dá)到應(yīng)有精度的目的。一、實(shí)驗(yàn)?zāi)康?、理解插值的基本概念，掌握各種插值方法，包括拉格朗日插值和牛頓插值等，注意其不同特點(diǎn)；2、通過實(shí)驗(yàn)進(jìn)一步理解并掌握各種插值的基本算法。二、實(shí)驗(yàn)任務(wù)1、 已知函數(shù)表 0.56160 0.56280 0.56401 0.56521 0.82741 0.82659 0.82577 0.82495用二次拉格朗日插值多項(xiàng)式求時的函數(shù)近似值。2、 已知函數(shù)表 0.4 0.55 0.65 0.8 0.9 0.41075 0.5781

5、5 0.69675 0.88811 1.02652用牛頓插值多項(xiàng)式求和。1.function y,R=lagranzi(X,Y,x,M)x=0.5635;M=2;X=0.56160,0.56280,0.56401,0.56521;Y=0.82741,0.82659,0.82577,0.82495;n=length(X); m=length(x);for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=abs(q1*(z-X(j); c1=c1*j

6、; end s=p*Y(k)+s; end y(i)=s;endR=M.*q1./c1;運(yùn)行結(jié)果：ans = 0.82612.N3(0.596)function y,R= newcz(X,Y,x,M) x=0.596; M=3; X=0.4,0.65,0.9; Y=0.41075,0.69675,1.02652; n=length(X); m=length(x); for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y' s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-

7、1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n); for k=(n-1):-1:1 C=conv(C,poly(X(k);d=length(C); C(d)=C(d)+A(k,k); end y(k)= polyval(C, z);endR=M*q1/c1;運(yùn)行結(jié)果：ans = 0.6313N4（0.895）function y,R= newcz(X,Y,x,M) x=0.895; M=4; X=0.4,0.55,0.65,0.8,0.9; Y=0.41075,0.57

8、815,0.69675,0.88811,1.02652; n=length(X); m=length(x); for t=1:m z=x(t); A=zeros(n,n);A(:,1)=Y' s=0.0; p=1.0; q1=1.0; c1=1.0; for j=2:n for i=j:n A(i,j)=(A(i,j-1)-A(i-1,j-1)/(X(i)-X(i-j+1); end q1=abs(q1*(z-X(j-1);c1=c1*j; end C=A(n,n);q1=abs(q1*(z-X(n); for k=(n-1):-1:1 C=conv(C,poly(X(k);d=len

9、gth(C); C(d)=C(d)+A(k,k); end y(k)= polyval(C, z);endR=M*q1/c1;運(yùn)行結(jié)果：ans = 1.0194實(shí)驗(yàn)三、解線性方程組的直接法解線性方程組的直接法是指經(jīng)過有限步運(yùn)算后能求得方程組精確解的方法。但由于實(shí)際計(jì)算中舍入誤差是客觀存在的，因而使用這類方法也只能得到近似解。目前較實(shí)用的直接法是古老的高斯消去法的變形，即主元素消去法及矩陣的三角分解法。引進(jìn)選主元的技巧是為了控制計(jì)算過程中舍入誤差的增長，減少舍入誤差的影響。一般說來，列主元消去法及列主元三角分解法是數(shù)值穩(wěn)定的算法，它具有精確度較高、計(jì)算量不大和算法組織容易等優(yōu)點(diǎn)，是目前計(jì)算機(jī)上解

10、中、小型稠密矩陣方程組可靠而有效的常用方法。一、實(shí)驗(yàn)?zāi)康?、 了解求線性方程組的直接法的有關(guān)理論和方法；2、 會編制列主元消去法、LU分解法的程序；3、 通過實(shí)際計(jì)算，進(jìn)一步了解各種方法的優(yōu)缺點(diǎn)，選擇合適的數(shù)值方法。 二、實(shí)驗(yàn)任務(wù)、 用列主元高斯消去法求解方程組.2、用矩陣直接三角分解法求解方程組，其中， .1.主程序：function RA,RB,n,X=liezhu(A,b)B=A b; n=length(b); RA=rank(A); ('求矩陣的秩')RB=rank(B);zhica=RB-RA;if zhica>0,disp('請注意：因?yàn)镽A=RB，所

11、以此方程組無解.')returnendif RA=RB if RA=ndisp('請注意：因?yàn)镽A=RB=n，所以此方程組有唯一解.') X=zeros(n,1); C=zeros(1,n+1); for p= 1:n-1Y,j=max(abs(B(p:n,p); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;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);endend b=B(1:n,n+1);A=B(1:n,1:n); X(n)=b(n)/A

12、(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('請注意：因?yàn)镽A=RB<n，所以此方程組有無窮多解.')endend計(jì)算程序：A=0.101 2.304 3.555;-1.347 3.712 4.623;-2.835 1.072 5.643; b=1.183;2.137;3.035; RA,RB,n,X=liezhu(A,b)運(yùn)行結(jié)果：ans =求矩陣的秩請注意：因?yàn)镽A=RB=n，所以此方程組有唯一解.RA = 3RB = 3n = 3X = -0.3982 0

13、.01380.33512. 程序：function X=LUjfcz(A,b)n,n =size(A); X=zeros(n,1); Y=zeros(n,1); C=zeros(1,n);r=1:n; for p=1:n-1max1,j=max(abs(A(p:n,p);C=A(p,:); A(p,:)= A(j+p-1,:); A(j+p-1,:)=C;g=r(p); r(p)= r(j+p-1);r(j+p-1)=g; if A(p,p)=0 disp('A是奇異陣,方程組無唯一解'); break;endfor k=p+1:nH= A(k,p)/A(p,p); A(k,p

14、) = H; A(k,p+1:n)=A(k,p+1:n)- H* A(p,p+1:n);endendY(1)=b(r(1);for k=2:nY(k)= b(r(k)- A(k,1:k-1)* Y(1:k-1); endX(n)= Y(n)/ A(n,n);for i=n-1:-1:1 X(i)= (Y(i)- A(i, i+1:n) * X (i+1:n)/ A(i,i);endEnd計(jì)算程序：A=1,2,-12,8;5,4,7,-2;-3,7,9,5;6,-12,-8,3; b=27;4;11;49;X=LUjfcz(A,b)運(yùn)行結(jié)果：X = 3.0000 -2.0000 1.00005.

15、0000實(shí)驗(yàn)四、解線性方程組的迭代法解線性方程組的迭代法是用某種極限過程去逐步逼近線性方程組精確解的方法，即是從一個初始向量出發(fā)，按照一定的迭代格式產(chǎn)生一個向量序列，使其收斂到方程組的解。迭代法的優(yōu)點(diǎn)是所需計(jì)算機(jī)存儲單元少，程序設(shè)計(jì)簡單，原始系數(shù)矩陣在計(jì)算過程中始終不變等。但迭代法存在收斂性及收斂速度問題。迭代法是解大型稀疏矩陣方程組的重要方法。一、實(shí)驗(yàn)?zāi)康?1、熟悉迭代法的有關(guān)理論和方法；2、會編制雅可比迭代法、高斯-塞德爾迭代法的程序；3、注意所用方法的收斂性及其收斂速度問題。二、實(shí)驗(yàn)任務(wù)1、用雅可比迭代法解方程組.注意：若用高斯-塞德爾迭代法則發(fā)散。2、用高斯-塞德爾迭代法解方程組.注意

16、：若用雅可比迭代法則發(fā)散。1.主程序：function X=jacdd(A,b,X0,P,wucha,max1)n m=size(A);for k=1:max1kfor j=1:mX(j)=(b(j)-A(j,1:j-1,j+1:m)*X0(1:j-1,j+1:m)/A(j,j);endXdjwcX=norm(X'-X0,P); xdwcX=djwcX/(norm(X',P)+eps); X0=X'if (djwcX<wucha)&(xdwcX<wucha)disp('請注意：雅可比迭代收斂,此方程組的精確解jX和近似解X如下：')r

17、eturnendendif (djwcX>wucha)&(xdwcX>wucha)disp('請注意：雅可比迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1 ')End計(jì)算程序：A=1 2 -2;1 1 1;2 2 1;b=7;2;5; X0=0 0 0'X=jacdd(A,b,X0,inf,0.01,100)運(yùn)行結(jié)果：k = 1X = 7 2 5k = 2X = 13 -10 -13k = 3X = 1 2 -1k = 4X = 1 2 -12.主程序：function X=gsdddy(A,b,X0,P,wucha,max1)D=diag(diag(A);U

18、=-triu(A,1);L=-tril(A,-1); dD=det(D);if dD=0disp('請注意：因?yàn)閷蔷仃嘍奇異，所以此方程組無解.')elsedisp('請注意：因?yàn)閷蔷仃嘍非奇異，所以此方程組有解.')iD=inv(D-L); B2=iD*U;f2=iD*b;jX=Ab; X=X0; n m=size(A);for k=1:max1X1= B2*X+f2; djwcX=norm(X1-X,P);xdwcX=djwcX/(norm(X,P)+eps);if (djwcX<wucha)|(xdwcX<wucha) return els

19、e k,X1',k=k+1;X=X1;endendif (djwcX<wucha)|(xdwcX<wucha) disp('請注意：高斯-塞德爾迭代收斂,此A的分解矩陣D,U,L和方程組的精確解jX和近似解X如下： ') elsedisp('請注意：高斯-塞德爾迭代的結(jié)果沒有達(dá)到給定的精度，并且迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1,方程組的精確解jX和迭代向量X如下： ')X=X'jX=jX'endendX=X'D,U,L,jX=jX'計(jì)算程序：A=1 0.9 0.9;0.9 1 0.9;0.9 0.9 1;b

20、=1.9;2.0;1.7;X0=0 0 0'X=gsdddy(A,b,X0,inf, 0.001,100)運(yùn)行結(jié)果：k =1ans =1.9000 0.2900 -0.2710k =2ans =1.8829 0.5493 -0.4890k = 3ans =1.8457 0.7789 -0.6622k =4ans =1.7949 0.9805 -0.7979k =5ans =1.7356 1.1560 -0.9025k =6ans =1.6718 1.3076 -0.9815k = 7ans =1.6065 1.4375 -1.0396k = 8ans =1.5419 1.5479 -1.0808k =9ans =1.4796 1.6411 -1.1086k =10ans =1.4208 1.7191 -1.1259k =11ans =1.3661 1.7838 -1.1349k =12