數值分析實驗報告

1、機電工程學院 機械工程 陳星星 6720150109 數值分析課程設計實驗報告實驗一 函數插值方法 一、問題提出 對于給定的一元函數的n+1個節(jié)點值。試用Lagrange公式求其插值多項式或分段二次Lagrange插值多項式。 數據如下： （1） 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.578150.696750.90 1.00 1.25382 求五次Lagrange多項式，計算,的值。（提示：結果為, ） 實驗步驟：第一步：先在matlab中定義lagran的M文件為拉格朗日函數代碼為：functionc,l=lagran(x,y)w=length(x)

2、;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;endc=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;>>p = lagran(x,y)回車得到：P = 121.6264 -422.7503 572.5667 -377.2549 1

3、21.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

4、.6264*x.5-422.7503*x.4+572.5667*x.3-377.2549*x.2+121.9718*x-15.084y =0.6257 得到f（0.596）=0.6257 同理得到f（0.99）=1.0542（2） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 試構造Lagrange多項式，和分段三次插值多項式，計算的,值。（提示：結果為, ）實驗步驟：第一步定義functionc,l=lagran(x,y)w=length(x);n=w-1;l=zeros(w,w);for k=1:n+1 v=1; for j

5、=1:n+1 if(k=j) v=conv(v,poly(x(j)/(x(k)-x(j); end end l(k,:)=v;endc=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;>> p=lagran(x,y)回車得到：由此得出所求拉格朗日多項式為p（x）=0.0001x6-0.0016x5+0.0186x4-0.1175x3+0.4419x2-0.9683x+0.9950第三步：在編輯窗口輸入如下命令：>

6、> 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=4304240283865561*x.6/73786976294838206464 - 7417128346304051*x.5/4611686018427387904 + 223*x.4/12000 - 2821*x.3/24000 + 994976512675275*x

7、.2/2251799813685248 - 19367*x./20000 + 199/200y =0.1648 得到f（1.8）=0.1648 同理得到f（6.15）=0.0013實驗結論： 插值是在離散數據的基礎上補插連續(xù)函數，使得這條連續(xù)曲線通過全部給定的離散數據點，它是離散函數逼近的重要方法，利用它可通過函數在有限個點處的取值狀況，估算出函數在其他點處的近似值。實驗二 函數逼近與曲線擬合 一、問題提出 從隨機的數據中找出其規(guī)律性，給出其近似表達式的問題，在生產實踐和科學實驗中大量存在，通常利用數據的最小二乘法求得擬合曲線。 在某冶煉過程中，根據統計數據的含碳量與時間關系，試求含碳量與時間

8、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 第一步先寫出線性最小二乘法的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);enda=f'*f;b=f'*y'c=ab;c=flipud(c);第二步在命令窗口輸入：>>lspoly(0,5,10,15,20,25

9、,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.20370.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í)行得到如下圖實驗結論： 分析復雜實驗數據時,常采用分段曲線擬合方法。

10、利用此方法在段內可以實現最佳逼近,但在段邊界上卻可能不滿足連續(xù)性和可導性。分段函數的光滑算法,給出了相應的誤差分析.給出了該方法在分段曲線擬合中的應用方法以及凸輪實驗數據自動分段擬合。實驗四 線方程組的直接解法一、問題提出 給出下列幾個不同類型的線性方程組，請用適當算法計算其解。 1、 設線性方程組 2、 設對稱正定陣系數陣線方程組 3、 三對角形線性方程組 實驗步驟：列主元高斯消去法的matlab的M文件程序function x,det,index=Gauss(A,b)% 求線形方程組的列主元Gauss消去法，其中，% A為方程組的系數矩陣；% b為方程組的右端項；% x為方程組的解；% d

11、et為系數矩陣A的行列式的值；% index為指標變量，index=0表示計算失敗，index=1表示計算成功。n,m=size(A); nb=length(b);% 當方程組行與列的維數不相等時，停止計算，并輸出出錯信息。if n=m error('The rows and columns of matrix A must be equal!'); return;end% 當方程組與右端項的維數不匹配時，停止計算，并輸出出錯信息if m=nb error('The columns of A must be equal the length of b!'); r

12、eturn;end% 開始計算，先賦初值index=1;det=1;x=zeros(n,1);for k=1:n-1 % 選主元 a_max=0; for i=k:n if abs(A(i,k)>a_max a_max=abs(A(i,k);r=i; end end if a_max<1e-10 index=0;return; end % 交換兩行 if r>k 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(

13、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);enddet=det*A(n,n);% 回代過程if abs(A(n,n)<1e-10 index=0;return;endfor k=n:-1:1 for j=k+1:n b(k)=b(k)-A(k,j)*x(j); end x(k)=b(k)/A(k,k);end然后在命令窗口輸入>> 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

14、 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.000

15、0 -1.00002.0000高斯-約當消去法maltab的M文件程序function x,flag=Gau_Jor(A,b)% 求線形方程組的列主元Gauss-約當法消去法，其中，% A為方程組的系數矩陣；% b為方程組的右端項；% x為方程組的解；n,m=size(A); nb=length(b);% 當方程組行與列的維數不相等時，停止計算，并輸出出錯信息。if n=m error('The rows and columns of matrix A must be equal!'); return;end% 當方程組與右端項的維數不匹配時，停止計算，并輸出出錯信息if m=

16、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 max1<1e-10 falg='failure'return; end % 交換兩行 if r>k for j=k:n z=A(k,j);A(k,j)=A

17、(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 endend% 輸出xfor 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

18、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實驗結論：用LU法，調用matlab中的函數lu中，L往往不是一個下三角，但可以直接計算不用它的結果來計算，不

19、用進行行變換。如果進行行變b也要變，這樣會很麻煩。實驗八 常微分方程初值問題數值解法一、基本題 科學計算中經常遇到微分方程（組）初值問題，需要利用Euler法，改進Euler法，Rung-Kutta方法求其數值解，諸如以下問題： （1） 分別取h=0.1,0.2,0.4時數值解。 初值問題的精確解。 （2） 用r=3的Adams顯式和預 - 校式求解 取步長h=0.1，用四階標準R-K方法求值。 （3） 用改進Euler法或四階標準R-K方法求解 取步長0.01,計算數值解，參考結果 。（4）利用四階標準R- K方法求二階方程初值問題的數值解 （I） (II) (III) (IV) 實驗步驟：function x,y=euler(fun,x0,xfinal,y0,n);if nargin<5,n=50;endh=(xfinal-x0)/n;x(1)=x0;y(1)=y0;for i=1:n; x(i+1)=x(i)+h; y(i+1)=y(i)+h*feval(fun,x(i),y(i);end實驗程序及分析 （） （1）、算法程序 function E =Euler_1(fun,x0,y0,xN,N) %