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

1、 山東建筑大學(xué)數(shù)值計(jì)算A實(shí)驗(yàn)報(bào)告二級(jí)學(xué)院： 理學(xué)院 專 業(yè)： 信息與計(jì)算科學(xué) 指導(dǎo)教師： XXX 班級(jí)學(xué)號(hào)： 信計(jì)XX 20081212XX 姓 名： XXXX 實(shí)驗(yàn)一 高斯消去法（1）題目：用選主元素法和高斯消去法求解下列方程組：已知方程增廣矩陣：a = 2 -1 0 0 6 -1 -3 -2 0 1 -1 3 -2 0 0 0 0 -3 5 1（2）原理： (1)高斯消去法：相對(duì)于約當(dāng)消去法而言，總的來說就是將增廣矩陣化為下三角陣。(2) 順序選主元素法:相對(duì)于高斯消去法的唯一不同是要先按當(dāng)前要排列元素所在列大小進(jìn)行排列。（3）設(shè)計(jì)思想： (1)高斯消去法：首先讓每一行的元素除以該行的主對(duì)

2、角線元素。然后利用此行使位于下一行主對(duì)角線以前的元素變?yōu)?，依次類推。 (2)選主元素法：在高斯消去法的基礎(chǔ)上，每次進(jìn)行化上三角陣之前，重新排列各方程的位置。（4）對(duì)應(yīng)程序：（a）高斯消去法function y=gauss1(a,b)%高斯順序消去法m,n=size(a);if m=n disp(輸入錯(cuò)誤，系數(shù)矩證陣只能是方陣)endif n=length(b) disp(輸入錯(cuò)誤,常數(shù)項(xiàng)的個(gè)數(shù)應(yīng)與方程的個(gè)數(shù)相同)endfor k=1:n-1 for i=k+1:n a(i,k)=a(i,k)/a(k,k); b(i)=b(i)-a(i,k)*b(k); for j=k+1:n if a(k,

3、k)=0 disp(主元素為零，消去法無法繼續(xù)) ; break; else a(i,j)=a(i,j)-a(i,k)*a(k,j); end end endendb(n)=b(n)/a(n,n);for i=(n-1):-1:1 w=0; for j=(i+1):n w=w+a(i,j)*b(j); end b(i)=(b(i)-w)/a(i,i); end y=b;（b）高斯列主元消去法function z=gauss2(a,b,ep) %高斯列主元素消元法if nargin=2 ep=0.000001endm,n=size(a);if m=n disp(輸入錯(cuò)誤，系數(shù)矩證陣只能是方陣)e

4、ndif n=length(b) disp(輸入錯(cuò)誤,常數(shù)項(xiàng)的個(gè)數(shù)應(yīng)與方程的個(gè)數(shù)相同)endfor k=1:n-1 p=a(k,k);I=k; for i=k:n if abs(a(i,k)abs(p) p=a(i,k);I=i; end end if pa = 2 -1 0 0 -1 -3 -2 0 -1 3 -2 0 0 0 -3 5 b = 6 1 0 1 gauss1(a,b)方程組的解為ans = 35/12 -1/6 -41/24 -33/40 （6）實(shí)驗(yàn)體會(huì)：主元消去法和高斯消去法的確是兩個(gè)非常鍛煉人編程的方法，在編寫程序時(shí)，需要使用的大量的循環(huán)和分支結(jié)構(gòu)，但無論是高斯消去法還是

5、高斯列主元法，它們的原理還算不難理解，通過變成能夠較好的理解它們。實(shí)驗(yàn)二 解線性方程組的迭代法2.1實(shí)驗(yàn)?zāi)康恼莆战饩€性方程組的雅可比迭代和高斯-塞德爾迭代算法；培養(yǎng)編程與上機(jī)調(diào)試能力.2.2實(shí)驗(yàn)要求： （1）選擇一種計(jì)算機(jī)語言（Matlab）設(shè)計(jì)出雅可比（Jacobi）Gauss-Seidel、SOR迭代法，迭代法的算法程序，并且選擇不同的迭代次數(shù)，觀察輸出結(jié)果； （2）利用Matlab求方程組2.3 實(shí)驗(yàn)內(nèi)容計(jì)算書上的習(xí)題（1）題目： P61例2.5.1a=20 2 3;1 8 1;2 -3 15; b=24;12;30; x0=0;0;0;（2）對(duì)應(yīng)程序：Jacobi迭代法：functio

6、n X=jacobi(a,b,X0,ep)%Jacobi迭代法求解方程組if nargin=3 ep=1.0e-6;elseif nargin=ep)&(kabs(p) p=a(i,k);I=i; end end if p=ep z=0; end if I=k for j=k:n w=a(k,j); a(k,j)=a(I,j);a(I,j)=w; end u=b(k);b(k)=b(I); b(I)=u; end for i=k+1:n a(i,k)=a(i,k)/a(k,k); b(i)=b(i)-a(i,k)*b(k); for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a

7、(k,j); end endendb(n)=b(n)/a(n,n);for i=(n-1):-1:1 w=0; for j=(i+1):n w=w+a(i,j)*b(j); end b(i)=(b(i)-w)/a(i,i);end % disp(方程組的解為);z=b;gauss2(a,b,ep)（3）實(shí)驗(yàn)結(jié)果： Jacobi迭代法：迭代次數(shù)為k = 6ans = 0.76736302083333 1.13841395833333 2.12535579166667高斯消元法：ans = 0.76735380732015 1.13840976020194 2.12536811106437（4）實(shí)

8、驗(yàn)體會(huì)：Jacobi迭代法和高斯消元法都能很好的解決方程組的求解問題，在上機(jī)程中遇到了也遇到了不少問題，但最后在老師的悉心輔導(dǎo)下都得到了很好的解答，這兩個(gè)程序使我明白了要變出好的程序就需要我們積極思考問題，勇于發(fā)現(xiàn)和解決問題。實(shí)驗(yàn)三 矩陣特征值問題計(jì)算3.1實(shí)驗(yàn)?zāi)康恼莆涨缶仃嚨奶卣髦岛椭魈卣飨蛄康膬绶?；培養(yǎng)編程與上機(jī)調(diào)試能力.3.2實(shí)驗(yàn)要求（1） 選擇一種計(jì)算機(jī)語言設(shè)計(jì)出冪法求主特征值和相應(yīng)特征向量的程序，并且選擇不同的初值，觀察所需的迭代次數(shù)和迭代結(jié)果.（2） 利用Matlab求特征值和特征向量 調(diào)用格式1： eig(A) %得到特征值列向量調(diào)用格式2：，其中為由特征列向量構(gòu)成的方陣，為由特

9、征值構(gòu)成的對(duì)角陣. %得到特征值和所對(duì)應(yīng)的特征向量3.3 實(shí)驗(yàn)內(nèi)容計(jì)算書上的習(xí)題（1）題目： P81例3.1.1A=2 4 6;3 9 15;4 16 36 x0=1;1;1（2）對(duì)應(yīng)程序： 乘冪法function y=chm(a,x0,k)%乘冪法求主特征值及特征向量if nargin=eps)&(kw w=abs(a(i,j); p=i; q=j; end endend if a(p,p)=a(q,q) theta=sign(a(p,q)*pi/4; G(p,p)=a(p,p)*(cos(theta)2+a(q,q)*(sin(theta)2+a(p,q)*sin(2*theta); G(

10、q,q)=a(p,p)*(sin(theta)2+a(q,q)*(cos(theta)2-a(p,q)*sin(2*theta); G(p,q)=(a(q,q)-a(p,p)*(cos(theta)*(sin(theta)+a(p,q)*(cos(2*theta); G(q,p)=G(p,q); for i=1:n R1(i,p)=R(i,p)*cos(theta)+R(i,q)*sin(theta); R1(i,q)=-R(i,p)*sin(theta)+R(i,q)*cos(theta); if (i=p)&(i=q) G(p,i)=a(i,p)*(cos(theta)+a(i,q)*(s

11、in(theta); G(i,p)=G(p,i); G(q,i)=-a(i,p)*(sin(theta)+a(i,q)*(cos(theta); G(i,q)=G(q,i); end end else gasi=(a(p,p)-a(q,q)/(2*a(p,q); t=sign(gasi)*(-abs(gasi)+sqrt(1+gasi2); G(p,p)=a(p,p)*(1/(1+t2)+a(q,q)*(t2/(1+t2)+a(p,q)*2*(t/(1+t2); G(q,q)=a(p,p)*(t2/(1+t2)+a(q,q)*(1/(1+t2)-a(p,q)*2*(t/(1+t2); G(p,

12、q)=(a(q,q)-a(p,p)*(t/(1+t2)+a(p,q)*(1/(1+t2)-(t2/(1+t2); G(q,p)=G(p,q); for i=1:n R1(i,p)=R(i,p)*(1/sqrt(1+t2)+R(i,q)*(t/sqrt(1+t2); R1(i,q)=-R(i,p)*(t/sqrt(1+t2)+R(i,q)*(1/sqrt(1+t2); if (i=p)&(i=q) G(p,i)=a(i,p)*(1/sqrt(1+t2)+a(i,q)*(t/sqrt(1+t2); G(i,p)=G(p,i); G(q,i)=-a(i,p)*(t/sqrt(1+t2)+a(i,q)

13、*(1/sqrt(1+t2); G(i,q)=G(q,i); end end end SS=0; for i=2:n for j=1:i-1 SS=SS+G(i,j)2; end end k=k+1;end disp(G的對(duì)角線即為近似的特征值) G disp(R的列向量為相應(yīng)的近似特征向量), RJeig(A,10e-9)（3）實(shí)驗(yàn)結(jié)果： 乘冪法：m = 43.87999781719956j = 3lanmda = 43.87999781719956u1 = 0.18586790191950 0.44603288809510 1.00000000000000經(jīng)典Jacobi法：G的對(duì)角線即為

14、近似的特征值G = 0.03838874652912 0.00000000000020 0.00000000013720 0.00000000000020 2.33485967315988 -0.00000086319675 0.00000000013720 -0.00000086319675 44.62675158031099R的列向量為相應(yīng)的近似特征向量R = 0.80459896338454 0.57043251684401 0.16500682363931 -0.58070234961462 0.69775944861021 0.41942404917606 0.12411804571

15、692 -0.43328800537539 0.89266803187144（4）實(shí)驗(yàn)體會(huì)：本次程序是由老師提供的，但是同學(xué)們都認(rèn)真閱讀過了，我發(fā)現(xiàn)此程序是相當(dāng)復(fù)雜的，要對(duì)矩陣進(jìn)行迭代等很多操作，我還了解到我們?cè)谔幚磔^大的問題是必須要使用到程序，因此讓我對(duì)程序產(chǎn)生的濃厚的興趣，同時(shí)我們也認(rèn)識(shí)到我們所編寫的程序穩(wěn)定性很差，因此我們還需要更多的練習(xí)。實(shí)驗(yàn)四 插值法4.1 實(shí)驗(yàn)?zāi)康?掌握插值法的基本思路和步驟； 培養(yǎng)編程與上機(jī)調(diào)試能力。4.2 實(shí)驗(yàn)要求用Matlab和插值中的某種具體算法編寫代碼并執(zhí)行，完成解決具體問題。4.3 實(shí)驗(yàn)內(nèi)容計(jì)算書上的習(xí)題Matlab SplineTools（1）題目：

16、x = 0 0.2000 0.4000 0.6000 0.8000y = 1.0000 1.2214 1.4918 1.8221 2.2255x0 = 0.1500（2）原理：（3）設(shè)計(jì)思想： （4）對(duì)應(yīng)程序：function s=Newton(x,y,x0,nn)%nn為Newton插值多項(xiàng)式的次數(shù)nx=length(x);ny=length(y);if nx=ny warning (向量x與y的長度應(yīng)該相同) return;endm=length(x0);for i=1:m yy=y; t=0; kk=1; while(kk=nn) kk=kk+1; for k=kk:nx yy(k)=(

17、yy(k)-yy(kk-1)/(x(k)-x(kk-1); end end t=yy(1); for k=2:nn+1 u=1.0; j=1; while(j Newton(x,y,x0,3)ans =1.1619 Newton(x,y,x0,4)ans =1.1618（7）實(shí)驗(yàn)體會(huì)：本實(shí)驗(yàn)通過對(duì)牛頓插值公式的程序化，是我們知道了可以用插值公式來解決函數(shù)問題，雖然我們只練習(xí)了牛頓公式，但是我們對(duì)插值公式有了很好的理解，對(duì)于在實(shí)驗(yàn)過程中遇到的問題都在老師的細(xì)心輔導(dǎo)下得到了很好的解決，本次試驗(yàn)受益匪淺。實(shí)驗(yàn)五 最小二乘法5.1 實(shí)驗(yàn)?zāi)康?掌握最小二乘法的基本思路和擬合步驟； 培養(yǎng)編程與上機(jī)調(diào)試能力

18、。5.2實(shí)驗(yàn)要求用Matlab和插值中的某種具體算法編寫代碼并執(zhí)行，完成解決具體問題。5.3 實(shí)驗(yàn)內(nèi)容計(jì)算書上的習(xí)題（1） 題目： 已知一組實(shí)驗(yàn)數(shù)據(jù)如下：X1234Y1.953.053.553.85求問題的最小二乘法。（2）原理：已知數(shù)據(jù)對(duì)，求多項(xiàng)式使得為最小，這就是一個(gè)最小二乘問題。（3）設(shè)計(jì)思想： 用線性函數(shù)為例，擬合給定數(shù)據(jù)。算法描述：步驟1：輸入值，及。步驟2：建立法方程組。步驟3：解法方程組。步驟4：輸出。（4）對(duì)應(yīng)程序：x=1:1:4;y=1.95,3.05,3.55,3.85;p=Polyfit(x,y,3);x1=1:0.1:4;y1=polyval(p,x1);plot(x,

19、y,*r,x1,y1,b)（5）實(shí)驗(yàn)結(jié)果： ans = 0.0667 -0.7000 2.7333 -0.1500（6）圖形（如果可視化）（7）實(shí)驗(yàn)體會(huì)：通過本次試驗(yàn)，我們發(fā)現(xiàn)MATLAB工具箱的功能是很強(qiáng)大的，MATLAB工具箱提供了最小二乘擬合的專門的函數(shù)，我們可以通過調(diào)用相應(yīng)的函數(shù)就可以達(dá)到我們需要的擬合結(jié)果，但要了解最小二乘擬合的思想，所以老師建議我們自己編寫相應(yīng)的代碼，通過本次試驗(yàn)，也極大的激發(fā)了我們的編程興趣，是我們認(rèn)識(shí)到了我們的不足之處，我們編寫的程序都是在正常的運(yùn)行情況下基本不出錯(cuò)，但是但出現(xiàn)異常時(shí)整個(gè)程序就會(huì)崩潰，也就是說我們編寫的程序不夠健壯，這需要我們?cè)谝院蟮膶W(xué)習(xí)中不斷進(jìn)

20、取。實(shí)驗(yàn)六 復(fù)化求積公式與龍貝格算法6.1 實(shí)驗(yàn)?zāi)康?掌握復(fù)化求積公式與龍貝格算法的基本思路和迭代步驟； 培養(yǎng)編程與上機(jī)調(diào)試能力。6.2 實(shí)驗(yàn)要求編寫程序，要求給出相應(yīng)習(xí)題的實(shí)驗(yàn)結(jié)果。6.3 實(shí)驗(yàn)內(nèi)容（1）題目： 用龍貝格算法計(jì)算：（2）原理：龍貝格算法利用外推法，提高了計(jì)算精度，加快了收斂速度。 對(duì)每一個(gè)從2做到，一直做到小于給定的精度是停止計(jì)算。其中（復(fù)化梯度求積公式），（3）設(shè)計(jì)思想： 步驟1：輸入?yún)^(qū)間端點(diǎn)，精度控制值，循環(huán)次數(shù)，定義函數(shù)，取， 步驟2：for to 步驟3：數(shù)據(jù)積分近似值。（4）對(duì)應(yīng)程序：function s=romberg1(fun,a,b,ep)% 龍貝格算法 -求

21、積公式% 此算法只外推到Romberg序列if nargin=3 ep=1.0e-6;elseif nargin=ep area=0.0; n=n+1; h=h/2; for i=1:m area=area+feval(fun,h*(2*i-1)+a); end t(n+1,1)=0.5*t(n,1)+area*h; m=2*m; if n4 for j=1:3 for i=1:n-j t(i,j+1)=(4(j)*t(i+1,j)-t(i,j)/(4(j)-1); end end t1=t(n-4,4); t2=t(n-3,4); endend disp(用Romberg序列求得積分近似值為

22、);s=t2;endfunction y=fun(x)if x=0 y=1;y=sin(x)/xromberg1(fun,0,1,10e-6)（5）實(shí)驗(yàn)結(jié)果： ans = 0.94608307036726（6）實(shí)驗(yàn)體會(huì)：以前不知道用程序怎么求積分，但通過本次實(shí)驗(yàn)我知道了利用程序是可以求積分的，它主要用的的有限次的迭代，將連續(xù)的無限多點(diǎn)的問題有線化，然后就可以利用計(jì)算機(jī)來解決積分問題，本程序由于比較復(fù)雜，所以老師給出了，雖然程序是由老師給出的，但是我們?cè)谡{(diào)用該程序時(shí)也基本明白了程序?qū)崿F(xiàn)的，在實(shí)驗(yàn)過程中我們仔細(xì)閱讀了程序，所以基本上可以再現(xiàn)，但我個(gè)人覺得這還遠(yuǎn)遠(yuǎn)不夠，我們需要自己獨(dú)立思考學(xué)習(xí)編寫，

23、不斷的提高自己的編程能力。實(shí)驗(yàn)七常微分方程數(shù)值解法7.1實(shí)驗(yàn)?zāi)康恼莆粘Ｎ⒎址匠虜?shù)值解的常用算法；培養(yǎng)編程與上機(jī)調(diào)試能力.求解常微分方程初值、邊值問題的解.7.2實(shí)驗(yàn)題目：下述方面的相應(yīng)習(xí)題（1）用改進(jìn)的歐拉公式，求解常微分方程初值問題的解 （2）用四階龍格-庫塔公式解初值問題（3）求解常微分方程邊值問題的解7.3實(shí)驗(yàn)要求（1） 選擇一種計(jì)算機(jī)語言設(shè)計(jì)出改進(jìn)歐拉法和四階龍格-庫塔法方法求解常微分方程初值問題的程序，觀察運(yùn)行結(jié)果.（2） 利用Matlab求解常微分方程初值問題 函數(shù)dsolve()用于求解微分方程.Dy表示：dy/dt（t 為缺省的自變量），Dny表示y對(duì)t的n階導(dǎo)數(shù). Matla

24、b6.5環(huán)境下操作如下： y=dsolve(Dy=y*y,y(0)=1) %求解題目1 y=dsolve(Dy=y/t,y(2.0)=1) %求解題目2 （3）ode23， ode45（1） 利用最小二乘法擬合通過改進(jìn)歐拉法求出微分方程的一系列數(shù)值解的近似函數(shù)方程.并利用Matlab的繪圖功能畫出函數(shù)的曲線。（1）題目1： y=x2+y2,0x1, y(0)=0, h=0.1題目2： y=3*y/(1+x),0x1, y(0)=1, h=0.2（2）原理：（3）設(shè)計(jì)思想： （4）對(duì)應(yīng)程序1：function s=Eulerc(fun,x0,xn,y0,h)%h為步長，改進(jìn)的求微分方程數(shù)值解的E

25、uler公式。if nargin5 error returnendn=(xn-x0)/h;for i=1:n t1=y0+h*feval(fun,x0,y0); x0=x0+h; t2=y0+h*feval(fun,x0,t1); y0=(t1+t2)/2; y(i)=y0;ends=y;return對(duì)應(yīng)程序2：function s=RK(fun,x0,xn,y0,h)%h為步長if nargin Eulerc(t2,x0,xn,y0,0.1)ans =1.1110 1.2515 1.4361 1.6880 2.0488 2.6000 3.5290 5.3715 10.3483 38.1343

26、實(shí)驗(yàn)結(jié)果2: RK(t3,x0,xn,y0,0.2)ans = 0 1.0000 0.2000 1.7275 0.4000 2.7430 0.6000 4.0942 0.8000 5.8292 1.0000 7.9960（6）圖形（如果可視化）（7）實(shí)驗(yàn)體會(huì)通過本次試驗(yàn)，我們領(lǐng)略到了程序的魅力，是我們對(duì)程序產(chǎn)生了濃厚的興趣，雖然在實(shí)驗(yàn)過程中我們遇到了比想象中要大的問題，但最后通過老師的悉心輔導(dǎo)和我們自己的不懈努力，最終都得到了圓滿的解答，本次實(shí)驗(yàn)讓我們看到了利用程序來求解微分方程的可能，我感覺到到程序也是一種藝術(shù)。實(shí)驗(yàn)八 非線性方程求根8.1實(shí)驗(yàn)?zāi)康模赫莆斩址?、牛頓迭代法等常用的非線性方程迭

27、代算法；培養(yǎng)編程與上機(jī)調(diào)試能力.8.2實(shí)驗(yàn)題目及參考結(jié)果：題目 求方程在1.5 附近的根.原方程的根為8.3實(shí)驗(yàn)要求：（1） 選擇一種計(jì)算機(jī)語言設(shè)計(jì)出二分法和牛頓法的程序，并且選擇不同的初值，觀察所需的迭代次數(shù)和迭代結(jié)果;（2） 分析二分法和牛頓法在非線性方程求根中的優(yōu)缺點(diǎn)和收斂速度二分法簡單易行，但只有線性收斂，且僅限于求實(shí)根；牛頓法也是一種簡單的迭代法，具有二階收斂速度（在單根鄰近處）的特點(diǎn)，但對(duì)初值的選擇比較苛刻，否則可能不收斂. (3) Matlab6.1環(huán)境下操作如下：f=1 1 -3 -3;roots(f)%多項(xiàng)式求根fx=x3+x2-3*x-3;fzero(fx,1.5)%非線性

28、方程式的實(shí)根（1）實(shí)驗(yàn)題目： 求方程在1.5 附近的根.原方程的根為（2）基本思想及原理：二分法對(duì)所給的非線性方程進(jìn)行初步確定其根所在區(qū)間為，此時(shí)有，令，判斷是否成立，成立則將，否則令，對(duì)此利用該方法，令逐步逼近方程的解。牛頓迭代法選取適當(dāng)?shù)某踔担?jì)算，然后過點(diǎn)作切線，則該切線與x軸的交點(diǎn)為，據(jù)此得到迭代公式為，當(dāng)時(shí)，得到方程的近似解。（3）對(duì)應(yīng)程序：二分法function c=midpart(a,b,eps,m) format longformat compactif myf(a)*myf(b)0 disp(該區(qū)間無根！) return;endwhile abs(a-b)=eps for k=1:m c=(a+b)/2 if myf(a)*myf(c)=eps for k=1:m x0=x1; x1=x0-dao(x0) k=k+1 end end disp(x1); function y=dao(x)y=(x3)+(x2)-(3*x)-3;y1=3*(x2)+2*x-3;y=y/y1;（4）實(shí)驗(yàn)結(jié)果： 二分法 midpart(1,2,0.00001,10)c = 1.50000000000000k = 2c = 1.75000000000000k = 3c = 1.62500000000000k =