《數(shù)值分析》課程實驗報告

《數(shù)值分析》課程實驗報告 精品文檔，僅供參考 《數(shù)值分析》課程實驗報告 《數(shù)值分析》課程實驗報告 姓 名： 學(xué) 號： 學(xué) 院： 機(jī) 電 學(xué) 院 日 期： 2015 年 X 月X 日 目 錄 實驗一 函數(shù)插值方法 1 實驗二 函數(shù)逼近與曲線擬合 5 實驗三 數(shù)值積分與數(shù)值微分 7 實驗四 線方程組的直接解法 9 實驗五 解線性方程組的迭代法 15 實驗六 非線性方程求根 19 實驗七 矩陣特征值問題計算 21 實驗八 常微分方程初值問題數(shù)值解法 24 實驗一 函數(shù)插值方法 一、問題提出 對于給定的一元函數(shù)的n+1個節(jié)點值。試用Lagrange公式求其插值多項式或分段二次Lagrange插值多項式。 數(shù)據(jù)如下： （1） 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多項式，和分段三次插值多項式，計算, 的值。（提示：結(jié)果為, ） （2） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 試構(gòu)造Lagrange多項式，計算的,值。（提示：結(jié)果為, ） 二、要求 1、利用Lagrange插值公式 編寫出插值多項式程序； 2、給出插值多項式或分段三次插值多項式的表達(dá)式； 3、根據(jù)節(jié)點選取原則，對問題（2）用三點插值或二點插值，其結(jié)果如何； 4、對此插值問題用Newton插值多項式其結(jié)果如何。Newton插值多項式如下： 其中： 三、目的和意義 1、學(xué)會常用的插值方法，求函數(shù)的近似表達(dá)式，以解決其它實際問題； 2、明確插值多項式和分段插值多項式各自的優(yōu)缺點； 3、熟悉插值方法的程序編制； 4、如果繪出插值函數(shù)的曲線，觀察其光滑性。 四、實驗步驟 （1） 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.57815 0.69675 0.90 1.00 1.25382 求五次Lagrange多項式，和分段三次插值多項式，計算, 的值。（提示：結(jié)果為, ） 第一步：先在matlab中定義lagran的M文件為拉格朗日函數(shù)代碼為： function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end 第二步：然后在matlab命令窗口輸入： >>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >> lagran(x,y) 回車得到： ans =121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845 由此得出所求拉格朗日多項式為 p（x）=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步：在編輯窗口輸入如下命令： >> x=[0.4 0.55 0.65 0.80,0.95 1.05]; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.0845; >> plot(x,y) 命令執(zhí)行后得到如下圖所示圖形，然后 >> x=0.596; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.6262 得到f（0.596）=0.6262 同理得到f（0.99）=1.0547 （2） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 試構(gòu)造Lagrange多項式，和分段三次插值多項式，計算的,值。（提示：結(jié)果為, ） 實驗步驟： 第一步定義 function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end 定義完拉格朗日M文件 第二步：然后在matlab命令窗口輸入： >>>> x=[1 2 3 4 5 6 7]; y=[0.368 0.135 0.050 0.018 0.007 0.002 0.001]; >> lagran(x,y) 回車得到： ans =0.0001 -0.0016 0.0186 -0.1175 0.4419 -0.9683 0.9950 由此得出所求拉格朗日多項式為 p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950 第三步：在編輯窗口輸入如下命令： >> x=[1 2 3 4 5 6 7]; >> y=0.0001*x.^6-0.0016*x.^5+0.0186*x.^4-0.1175*x.^3+0.4419*x.^2-0.9683*x+0.9950; >> plot(x,y) 命令執(zhí)行后得到如下圖所示圖形，然后 >> x=1.8; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718*x-15.084 y =0.1650 得到f（0.596）=0.6262 同理得到f（6.15）=2.3644 五、實驗結(jié)論 插值是在離散數(shù)據(jù)的基礎(chǔ)上補(bǔ)插連續(xù)函數(shù)，使得這條連續(xù)曲線通過全部給定的離散數(shù)據(jù)點，它是離散函數(shù)逼近的重要方法，利用它可通過函數(shù)在有限個點處的取值狀況，估算出函數(shù)在其他點處的近似值。 實驗二 函數(shù)逼近與曲線擬合 一、問題提出 從隨機(jī)的數(shù)據(jù)中找出其規(guī)律性，給出其近似表達(dá)式的問題，在生產(chǎn)實踐和科學(xué)實驗中大量存在，通常利用數(shù)據(jù)的最小二乘法求得擬合曲線。 在某冶煉過程中，根據(jù)統(tǒng)計數(shù)據(jù)的含碳量與時間關(guān)系，試求含碳量與時間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 二、要求 1、用最小二乘法進(jìn)行曲線擬合； 2、近似解析表達(dá)式為； 3、打印出擬合函數(shù)，并打印出與的誤差，； 4、另外選取一個近似表達(dá)式，嘗試擬合效果的比較； 5、* 繪制出曲線擬合圖。 三、目的和意義 1、掌握曲線擬合的最小二乘法； 2、最小二乘法亦可用于解超定線代數(shù)方程組； 3、探索擬合函數(shù)的選擇與擬合精度間的關(guān)系 四、實驗步驟： 第一步先寫出線性最小二乘法的M文件 function c=lspoly(x,y,m) n=length(x); b=zeros(1:m+1); f=zeros(n,m+1); for k=1:m+1 f(:,k)=x.^(k-1); end a=f＇*f; b=f＇*y＇; c=a＼b; c=flipud(c); 第二步在命令窗口輸入： >>lspoly([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],2) 回車得到： ans = -0.0024 0.2037 0.2305 即所求的擬合曲線為y=-0.0024x2+0.2037x+0.2305 在編輯窗口輸入如下命令： >> x=[0,5,10,15,20,25,30,35,40,45,50,55]; >> y=-0.0024*x.^2+0.2037*x+0.2305; >> plot(x,y) 命令執(zhí)行得到如下圖 五、實驗結(jié)論 分析復(fù)雜實驗數(shù)據(jù)時,常采用分段曲線擬合方法。利用此方法在段內(nèi)可以實現(xiàn)最佳逼近,但在段邊界上卻可能不滿足連續(xù)性和可導(dǎo)性。分段函數(shù)的光滑算法,給出了相應(yīng)的誤差分析.給出了該方法在分段曲線擬合中的應(yīng)用方法以及凸輪實驗數(shù)據(jù)自動分段擬合。 實驗三 數(shù)值積分與數(shù)值微分 一、問題提出 選用復(fù)合梯形公式，復(fù)合Simpson公式，Romberg算法，計算 （1） （2） （3） （4） 二、要求 1、編制數(shù)值積分算法的程序； 2、分別用兩種算法計算同一個積分，并比較其結(jié)果； 3、分別取不同步長，試比較計算結(jié)果（如n = 10, 20等）； 4、給定精度要求ε，試用變步長算法，確定最佳步長。 三、目的和意義 1、深刻認(rèn)識數(shù)值積分法的意義； 2、明確數(shù)值積分精度與步長的關(guān)系； 3、根據(jù)定積分的計算方法，可以考慮二重積分的計算問題。 四、實驗步驟 第一步：編寫各種積分的程序 復(fù)合梯形程序如下： function I=TX(x,y) n=length(x);m=length(y); if n~=m error(＇The lengths of X and Y must be equal＇); return; end h=(x(n)-x(1))/(n-1); a=[1 2*ones(1,n-2) 1]; I=h/2*sum(a.*y); 復(fù)合Simpson程序如下： function s = simpr1(f,a,b,n) h=(b-a)/(2*n); s1=0; s2=0; for k=1:10 x=a+h*(2*k-1); s1=s1+feval(f,x); end for k=1:(10-1) x=a+h*2*k; s2=s2+feval(f,x); end s=h*(feval(f,a)+feval(f,b)+4*s1+2*s2)/3; end Romberg程序如下： function I = Romber_yang(fun,a,b,ep) if nargin1; T(m+1,k+1)=(4^k*T(m+1,k)-T(m,k))/(4^k-1); M=M/2;k=k+1; end if abs(T(k,k)-T(k-1,k-1))a_max a_max=abs(A(i,k));r=i; end end if a_maxk for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end % 消元過程 for i=k+1:n m=A(i,k)/A(k,k); for j=k+1:n A(i,j)=A(i,j)-m*A(k,j); end b(i)=b(i)-m*b(k); end det=det*A(k,k); end det=det*A(n,n); % 回代過程 if abs(A(n,n))> A=[4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5 7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1]; >> b=[5 12 3 2 3 46 13 38 19 -21]; >> gauss(A,b) ans = 1.0000 -1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000 -1.0000 2.0000 高斯-約當(dāng)消去法maltab的M文件程序 function [x,flag]=Gau_Jor(A,b) % 求線形方程組的列主元Gauss-約當(dāng)法消去法，其中， % A為方程組的系數(shù)矩陣； % b為方程組的右端項； % x為方程組的解； [n,m]=size(A); nb=length(b); % 當(dāng)方程組行與列的維數(shù)不相等時，停止計算，并輸出出錯信息。 if n~=m error(＇The rows and columns of matrix A must be equal!＇); return; end % 當(dāng)方程組與右端項的維數(shù)不匹配時，停止計算，并輸出出錯信息 if m~=nb error(＇The columns of A must be equal the length of b!＇); return; end % 開始計算，先賦初值 flag=＇ok＇;x=zeros(n,1); for k=1:n % 選主元 max1=0; for i=k:n if abs(A(i,k))>max1 max1=abs(A(i,k));r=i; end end if max1k for j=k:n z=A(k,j);A(k,j)=A(r,j);A(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z; end % 消元過程 b(k)=b(k)/A(k,k); for j=k+1:n A(k,j)=A(k,j)/A(k,k); end for i=1:n if i~=k for j=k+1:n A(i,j)=A(i,j)-A(i,k)*A(k,j); end b(i)=b(i)-A(i,k)*b(k); end end end % 輸出x for i=1:n x(i)=b(i); end 然后保存后在命令窗口輸入： >> A=[4 2 -4 0 2 4 0 0;2 2 -1 -2 1 3 2 0;-4 -1 14 1 -8 -3 5 6;0 -2 1 6 -1 -4 -3 3;2 1 -8 -1 22 4 -10 -3;4 3 -3 -4 4 11 1 -4;0 2 5 -3 -10 1 14 2;0 0 6 3 -3 -4 2 19]; >> b=[0 -6 20 23 9 -22 -15 45]; >> Gau_Jor(A,b) ans = 121.1481 -140.1127 29.7515 -60.1528 10.9120 -26.7963 5.4259 -2.0185 五、實驗結(jié)論 用LU法，調(diào)用matlab中的函數(shù)lu中，L往往不是一個下三角，但可以直接計算不用它的結(jié)果來計算，不用進(jìn)行行變換。如果進(jìn)行行變b也要變，這樣會很麻煩。 實驗五 解線性方程組的迭代法 一、問題提出 對實驗四所列目的和意義的線性方程組，試分別選用Jacobi 迭代法，Gauss-Seidel迭代法和SOR方法計算其解。 二、要求 1、體會迭代法求解線性方程組，并能與消去法做以比較； 2、分別對不同精度要求，如由迭代次數(shù)體會該迭代法的收斂快慢； 3、對方程組2，3使用SOR方法時，選取松弛因子ω=0.8，0.9，1，1.1，1.2等，試看對算法收斂性的影響，并能找出你所選用的松弛因子的最佳者； 4、給出各種算法的設(shè)計程序和計算結(jié)果。 三、目的和意義 1、通過上機(jī)計算體會迭代法求解線性方程組的特點，并能和消去法比較； 2、運用所學(xué)的迭代法算法，解決各類線性方程組，編出算法程序； 3、體會上機(jī)計算時，終止步驟或k >（給予的迭代次數(shù)），對迭代法斂散性的意義； 4、體會初始解，松弛因子的選取，對計算結(jié)果的影響。 四、實驗步驟 第一步編寫實驗所需的Jacobi迭代法，Gauss-Seidel迭代法，SOR迭代法的程序。 Jacobi迭代法： function [x,k,index]=J(A,b,ep,itmax) if nargin Euler_1(＇f＇,0,3,2,20) 結(jié)果 ： >> Euler_1(＇f＇,0,3,2,20) T = 0 3.0000 0.1000 2.9836 0.2000 2.9517 0.3000 2.9058 0.4000 2.8481 0.5000 2.7810 0.6000 2.7073 0.7000 2.6297 0.8000 2.5511 0.9000 2.4739 1.0000 2.4004 1.1000 2.3325 1.2000 2.2714 1.3000 2.2177 1.4000 2.1717 1.5000 2.1332 1.6000 2.1017 1.7000 2.0765 1.8000 2.0567 1.9000 2.0414 2.0000 2.0299 五、實驗結(jié)論 很多科學(xué)技術(shù)和工程問題常用微分方程的形式建立數(shù)學(xué)模型，因此微分方程的求解是很有意義的，但對于絕大多數(shù)的微分方程問題很難或者根本不可能得到它的解析解，因此，研究微分方程的數(shù)值方法是非常有意義的。