1、精選優質文檔-傾情為你奉上數值分析實驗報告學院：電氣工程與自動化學院專業：控制理論與控制工程姓名：李亞學號：2014 年 12 月24日實驗一 函數插值方法 一、目的和意義 1、 學會常用的插值方法，求函數的近似表達式，以解決其它實際問題； 2、 明確插值多項式和分段插值多項式各自的優缺點； 3、 熟悉插值方法的程序編制； 4、 如果繪出插值函數的曲線，觀察其光滑性。二、實驗原理 1、 Lagrange插值公式 編寫出插值多項式程序； 2、 給出插值多項式或分段三次插值多項式的表達式； 三、實驗要求對于給定的一元函數的n+1個節點值。試用Lagrange公式求其插值多項式或分段二次Lagran

2、ge插值多項式。數據如下： （1） 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.578150.696750.90 1.00 1.25382 求五次Lagrange多項式，計算,的值。（提示：結果為, ） （2） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 試構造Lagrange多項式，和分段三次插值多項式，計算的,值。（提示：結果為, ）四、實驗過程1進入matlab開發環境；2根據實驗內容和要求編寫程序，程序如下所示，程序通過運用function函數編寫，生成.m文件。調用時只需要在命令窗口

3、調用y=Lagrange(A,input)就可以實現任意次數拉格朗日插值法求解。function y=Lagrange(A,input)a,b=size(A);x=input;y=0;for j=1:a Mj=1; Nj=1; for k=1:a if(k=j) continue; end Mj=Mj*(x-A(k,1); Nj=Nj*(A(j,1)-A(k,1); end y=y+A(j,2)*Mj/Nj;end3調試程序并運行程序；調用拉格朗日腳本文件對以上兩個表格數據求解，表格一對應MATLAB向量A;表格二對應向量I。在命令窗口調用y=Lagrange(A,input)，求解如下面截圖

4、。圖1 表一數據的解圖2 表二數據的解4實驗總結 通過對插值法算法編程，加深了對插值方法的理解，熟悉了MATLAB編寫腳本函數。通過計算機求解，能更加方便快捷求解。 實驗二 函數逼近與曲線擬合 一、目的和意義 1、掌握曲線擬合的最小二乘法； 2、最小二乘法亦可用于解超定線代數方程組； 3、探索擬合函數的選擇與擬合精度間的關系。二、實驗原理 對于給定的測量數據(xi,fi)(i=1,2,，n)，設函數分布為特別的，取為多項式 (j=0, 1,，m)則根據最小二乘法原理，可以構造泛函令 (k=0, 1,，m)則可以得到法方程求該解方程組，則可以得到解，因此可得到數據的最小二乘解三、實驗要求 1、用

5、最小二乘法進行曲線擬合； 2、近似解析表達式為；3、打印出擬合函數，并打印出與的誤差，； 4、另外選取一個近似表達式，嘗試擬合效果的比較； 5、繪制出曲線擬合圖。 四、實驗步驟： 1進入matlab開發環境；2根據實驗內容和要求編寫程序如下；代碼一公式S(x)=a1*t+a2*t2+a3*t3；代碼二公式S(x)=a2*t2+a3*t3+a4*t4.代碼一：function error=mintwomultiply(A)%S(x)=a1*t+a2*t2+a3*t3a,b=size(A);M=zeros(3);N=zeros(3,1);error=0;for i=1:a M(1,1)=M(1,1

6、)+A(i,1)*A(i,1); M(1,2)=M(1,2)+A(i,1)*A(i,1)2; M(2,1)=M(1,2); M(1,3)=M(1,3)+A(i,1)*A(i,1)3; M(3,1)=M(1,3); M(2,2)=M(2,2)+A(i,1)2*A(i,1)2; M(2,3)=M(2,3)+A(i,1)2*A(i,1)3; M(3,2)=M(2,3); M(3,3)=M(3,3)+A(i,1)3*A(i,1)3; N(1,1)=N(1,1)+A(i,1)*A(i,2); N(2,1)=N(2,1)+A(i,1)2*A(i,2); N(3,1)=N(3,1)+A(i,1)3*A(i,

7、2);end%a1,a2,a3=solve(M,N)I=MN;for i=1:aA(i,3)=I(1,1)*A(i,1)+I(2,1)*A(i,1)2+I(3,1)*A(i,1)3; error=error+(A(i,3)-A(i,2)2;endhold on;plot(A(:,1),A(:,3),'r','LineWidth',2);plot(A(:,1),A(:,2),'b','LineWidth',2);legend('原始圖像',擬合圖像',2);hold off;代碼二：function erro

8、r=mintwomultiply2(A)%S(x)=a2*t2+a3*t3+a4*t4a,b=size(A);M=zeros(3);N=zeros(3,1);error=0;for i=1:a M(1,1)=M(1,1)+A(i,1)2*A(i,1)2; M(1,2)=M(1,2)+A(i,1)2*A(i,1)3; M(2,1)=M(1,2); M(1,3)=M(1,3)+A(i,1)2*A(i,1)4; M(3,1)=M(1,3); M(2,2)=M(2,2)+A(i,1)3*A(i,1)3; M(2,3)=M(2,3)+A(i,1)3*A(i,1)4; M(3,2)=M(2,3); M(3

9、,3)=M(3,3)+A(i,1)4*A(i,1)4; N(1,1)=N(1,1)+A(i,1)2*A(i,2); N(2,1)=N(2,1)+A(i,1)3*A(i,2); N(3,1)=N(3,1)+A(i,1)4*A(i,2);end%a1,a2,a3=solve(M,N)I=MN;for i=1:a A(i,3)=I(1,1)*A(i,1)2+I(2,1)*A(i,1)3+I(3,1)*A(i,1)4 error=error+(A(i,3)-A(i,2)2;endhold on;plot(A(:,1),A(:,3),'r','LineWidth',2);

10、plot(A(:,1),A(:,2),'b','LineWidth',2);legend('原始圖像',擬合圖像',2);hold off;3調試程序并運行程序；在某冶煉過程中，根據統計數據的含碳量與時間關系，試求含碳量與時間t的擬合曲線。t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 根據以上數據調用代碼一和代碼二，結果如下圖，返回的結果為error為誤差，擬合圖像藍的為原始圖像，紅的為擬合圖像，

11、具體圖像見下圖:圖1 公式一誤差圖2 公式一擬合圖圖3 公式二誤差圖4 公式二擬合圖4實驗總結從不同公式產生的誤差來看，公式一明顯誤差小；從產生的擬合圖像來看公式一更穩定，更平穩。由此可得，公式的選取對結果的影響很大。實驗三 線方程組的直接解法一、目的和意義 1、通過該課題的實驗，體會模塊化結構程序設計方法的優點； 2、運用所學的計算方法，解決各類線性方程組的直接算法； 3、提高分析和解決問題的能力，做到學以致用；4、學習Gauss消去法的原理；5、了解列主元的意義；6、確定什么時候系數陣要選主元。二 實驗數學原理：由于一般線性方程在使用Gauss消去法求解時，從求解的過程中可以看到，若=0，

12、則必須進行行交換，才能使消去過程進行下去。有的時候即使0，但是其絕對值非常小，由于機器舍入誤差的影響，消去過程也會出現不穩定得現象，導致結果不正確。因此有必要進行列主元技術，以最大可能的消除這種現象。這一技術要尋找行r，使得并將第r行和第k行的元素進行交換，以使得當前的的數值比0要大的多。這種列主元的消去法的主要步驟如下：1. 消元過程對k=1,2,n-1,進行如下步驟。1) 選主元，記若很小，這說明方程的系數矩陣嚴重病態，給出警告，提示結果可能不對。2) 交換增廣陣A的r，k兩行的元素。 (j=k,n+1)3) 計算消元 (i=k+1,n; j=k+1,n+1)2. 回代過程對k= n, n

13、-1,1,進行如下計算至此，完成了整個方程組的求解。三、實驗要求 1、對下述方程組分別利用Gauss順序消去法與Gauss列主元消去法；2、應用結構程序設計編出通用程序； 3、比較計算結果，分析數值解誤差的原因； 4、設線性方程組如下 方程的解：四、實驗要求 1進入matlab開發環境；2根據實驗內容和要求編寫程序；代碼分為兩個，Gauss順序消去法與Gauss列主元消去法；代碼一：Gauss順序消去法function x=orderGauss(A,B)m,n=size(A);k=rank(A);if m=n&&m=k for i=1:n-1 for j=i+1:n B(j,:

14、)=B(j,1)-A(j,i)*B(i,1)/A(i,i); A(j,:)=A(j,:)-A(j,i)*A(i,:)/A(i,i); end end x(n,1)=B(n,1)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*x(j,1); end x(i,1)=(B(i,1)-sum)/A(i,i); endend代碼二：Gauss列主元消去法此代碼分為兩部分，即兩個功能函數。rows,columns,temp=findmax(A)和x=Gauss(A,B)。前面函數尋找向量中的最大值，返回所在的行，列，具體值。后面再find

15、max基礎進行列主元迭代，代碼如下：第一部分：尋找向量，矩陣最大值function rows,columns,temp=findmax(A)m n=size(A);if m=1&&n=1 rows=1; columns=1; temp=A(1,1); for i=2:m if abs(A(i,1)>abs(temp) temp=A(i,1); rows=i; end endelseif m=1&&n=1 rows=1; columns=1; temp=A(1,1); for i=2:n if abs(A(1,i)>abs(temp) temp=A(1

16、,i); columns=i; end end else rows=1; columns=1; temp=A(1,1); for i=1:m for j=1:n if abs(A(i,j)>abs(temp) temp=A(i,j); rows=i; columns=j; end end endend第二部分：高斯消元法求解function x=Gauss(A,B)m,n=size(A);k=rank(A);if m=n&&m=k for i=1:n-1 a,c=findmax(A(i:m,i); a=a+i-1; if c=0 continue; end if a=i

17、temp=A(i,:); A(i,:)=A(a,:); A(a,:)=temp; temp1=B(i,1); B(i,1)=B(a,1); B(a,1)=temp1; end for j=i+1:n B(j,:)=B(j,1)-A(j,i)*B(i,1)/c; A(j,:)=A(j,:)-A(j,i)*A(i,:)/c; end end x(n,1)=B(n,1)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*x(j,1); end x(i,1)=(B(i,1)-sum)/A(i,i); endend3調試程序并運行程序；用G

18、auss順序消去法與Gauss列主元消去法求解方程組的解，結果如下圖： 圖1 Gauss順序消去法的解 圖2 Gauss列主元消去法的解4實驗總結對比試驗結果，兩種方法的解和方程的解都一直。這是由于計算機有足夠的精度，如果沒有足夠的精度，誤差就容易被放大。產生不正確的結果。列主元消去法可以克服這種缺陷，求出相對精確的解。實驗四 解線性方程組的迭代法一、目的和意義 1、通過上機計算體會迭代法求解線性方程組的特點，并能和消去法比較； 2、運用所學的迭代法算法，解決各類線性方程組，編出算法程序； 3、體會上機計算時，終止步驟或k >（給予的迭代次數），對迭代法斂散性的意義； 4、體會初始解，松

19、弛因子的選取，對計算結果的影響。二、實驗原理設線性方程組 (1)的系數矩陣A可逆且主對角元素均不為零,令 并將A分解成 (2)從而(1)可寫成 令 其中. (3)以為迭代矩陣的迭代法(公式) (4)稱為雅可比(Jacobi)迭代法(公式),用向量的分量來表示,(4)為 (5)其中為初始向量.三、要求 1、體會迭代法求解線性方程組，并能與消去法做以比較；四、實驗過程 1進入matlab開發環境2根據實驗內容和要求編寫程序此代碼分為兩部分，即兩個功能函數。y=Jacobicondition(A)和X=Jacobi(A,B,X,input)。前面函數為判斷條件，判斷矩陣A是否可以進行雅克比迭代。后面X=Jacobi(A,B,X,input)進行迭代求解，具體代碼如下：代碼一：function y=Jacobicondition