張揚 2010212420數(shù)值計算方法實驗

1、 實驗一 插值法一、 算法設計思想:在插值計算中，取樣點的多少往往會影響所得插值函數(shù)優(yōu)化程度。一般情況下，取樣點越多所得插值函數(shù)越優(yōu)化，對應的函數(shù)值與標準函數(shù)值越接近。在MATLAB中實現(xiàn)分段線性插值，分段二次插值，拉格朗日插值的命令為interp1,其調(diào)用格式為: Y1=interp1(X,Y,X1，method)函數(shù)根據(jù)X,Y的值，計算函數(shù)在X1處的值。X,Y是兩個等長的已知向量，分別描述采樣點和樣本值，X1是一個向量或標量，描述欲插值點，Y1是一個與X1等長的插值結果。二、 程序清單：1. 分段線性插值#include<stdio.h>#include<math.h&g

2、t; double Lagrange2(double *x, double *y, double input) double output; int i; for (i=0;i<5;i+) if (xi <= input && xi+1 >= input) output=yi +(yi+1-yi)*(input-xi)/(xi+1-xi); break; return output;2.分段二次插值double Lagrange3(double *x,double *y,double u) int i,k=0; double v; for(i=0;i<6

3、;i+) if(u<x1) k=0; v=yk*(u-xk+1)*(u-xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); if(xi<u&&u<=xi+1)&&(fabs(u-xi)<=fabs(u-xi+1) k=i-1; v=yk*(u-xk+1)*(u-xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-x

4、k)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); if (xi<u&&u<=xi+1)&&fabs(u-xi)>fabs(u-xi+1) k=i; v=yk*(u-xk+1)*(u-xk+2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); if(u>x4) k=3; v=yk*(u-xk+1)*(u-xk+

5、2)/(xk-xk+1)*(xk-xk+2)+yk+1*(u-xk)*(u-xk+2)/(xk+1-xk)*(xk+1-xk+2)+yk+2*(u-xk)*(u-xk+1)/(xk+2-xk)*(xk+2-xk+1); return v;void main() double x6 = 0.0, 0.1, 0.195, 0.3, 0.401, 0.5,y6 = 0.39894,0.39695,0.39142,0.38138,0.36812,0.35206; double u; scanf("%lf",&u);拉格朗日插值：double Lagrange1(double

6、 *x, double *y, double xx) int i,j; double *a,yy=0.000; a=new double6; for(i=0;i< 6;i+) ai=yi; for(j=0;j< 6;j+) if(j!=i) ai*=(xx-xj)/(xi-xj); yy+=ai; delete a; return yy;double Lagrange2(double *x, double *y, double input) double output; int i; for (i=0;i<5;i+) if (xi <= input &&

7、 xi+1 >= input) output=yi +(yi+1-yi)*(input-xi)/(xi+1-xi); break; return output;三、 運行結果：四、 四種插值的比較：從以上對比函數(shù)圖象可以看出，分段線性插值其總體光滑程度不夠。在數(shù)學上，光滑程度的定量描述是函數(shù)(曲線) 的k階導數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性。一般情況下，階數(shù)越高光滑程度越好。分段線性插值具有零階光滑性，也就是不光滑。二次插值就是較低次數(shù)的多項式而達到較高階光滑性的方法。總體上分段線性插值具有以下特點：五、經(jīng)驗總結：分段線性插值在計算上具有簡潔方便的特點。分段線性插值與二次插值函數(shù)在

8、每個小區(qū)間上相對于原函數(shù)都有很強的收斂性，（舍入誤差影響不大），數(shù)值穩(wěn)定性好且容易在計算機上編程實現(xiàn)等優(yōu)點，但分段線性插值在節(jié)點處具有不光滑性的缺點(不能保證節(jié)點處插值函數(shù)的導數(shù)連續(xù))，從而不能滿足某些工程技術上的要求。實驗二 曲線擬合的最小二乘法一、 算法設計思想：對給定數(shù)據(jù) (i=0 ,1,2,3,.,m),一共m+1個數(shù)據(jù)點，取多項式P(x),使二、函數(shù)P(x)稱為擬合函數(shù)或最小二乘解，令似的 使得 其中，a0，a1,a2,an為待求未知數(shù)，n為多項式的最高次冪，由此，該問題化為求的極值問題。由多元函數(shù)求極值的必要條件： j=0,1,n得到： j=0,1,n這是一個關于a0，a1,a2,

9、an的線性方程組，用矩陣表示如下：因此，只要給出數(shù)據(jù)點 及其個數(shù)m，再給出所要擬合的參數(shù)n，則即可求出未知數(shù)矩陣（a0，a1,a2,an）二、程序清單： x=-1.0:0.5:2.0;y=-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552;plot(x,y,'*')xlabel 'x軸'ylabel 'y軸'title '散點圖'hold onm=6;n=3;A=zeros(n+1);for j=1:n+1 for i=1:n+1 for k=1:m+1 A(j,i)=A(j,i)+x(k)(

10、j+i-2) end endend;B=0 0 0 0;for j=1:n+1 for i=1:m+1 B(j)=B(j)+y(i)*x(i)(j-1) endendB=B'a=inv(A)*B;x=-1.0:0.0001:2.0;z=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3;plot(x,z) legend('離散點','y=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3')title('擬合圖')三、運行的結果：五、經(jīng)驗總結：在運籌學、統(tǒng)計學、逼近論和控制論中，最小二乘法都是很重要的求解方法。例如，它是統(tǒng)

11、計學中估計回歸參數(shù)的最基本方法。在實際問題中，怎樣由測量的數(shù)據(jù)設計和確定“最貼近”的擬合曲線？關鍵在選擇適當?shù)臄M合曲線類型，有時根據(jù)專業(yè)知識和工作經(jīng)驗即可確定擬合曲線類型；在對擬合曲線一無所知的情況下，不妨先繪制數(shù)據(jù)的粗略圖形，或許從中觀測出擬合曲線的類型；更一般地，對數(shù)據(jù)進行多種曲線類型的擬合，并計算均方誤差，用數(shù)學實驗的方法找出在最小二乘法意義下的誤差最小的擬合函數(shù)。實驗三 矩陣的特征值和特征向量一、算法設計思想：設是階方陣，若存在數(shù)和維非零向量使關系式 （1）成立，則稱這樣的數(shù)為方陣A的特征值；非零向量稱為A的對應于特征值的特征向量將（1）式改寫成， （2）得到一個含個未知數(shù)個方程的齊次

12、線性方程組，它有非零解的充分必要條件是其系數(shù)行列式一、即這個以為未知數(shù)的一元次方程稱為方陣的特征方程，其左端是的次多項式，記作，稱為方陣的特征多項式，顯然，的特征值就是特征方程的根；在復數(shù)范圍內(nèi)，階方陣有個特征值。從定義1不難看出，若非零列向量是方陣的對應于特征值的特征向量，則（為常數(shù)）也是的對應于的特征向量，即對應于一個特征值有無窮多個特征向量。反之，不同的特征值所對應的特征向量不相等，即一個特征向量只能屬于一個特征值。由上面的討論得到方陣的特征值與特征向量的求法：1）寫出方陣的特征方程，它的根就是的全部特征值。2）對每個特征值，齊次線性方程組（2）的每一個非零解都是的對應于的特征向量，只要

13、求出（2）的一個基礎解系，它們的線性組合（為非零向量時）就是的對應于的全部的特征向量。二、 程序清單：clc;clear;close; A=4,1,-1;16,-2,-2;2,-3,-1; X,B=eig(A)三、運行結果：A = 3 -1 -2 2 0 -2 2 -1 -1>> V,D=eig(A)V = 0.7276 -0.5774 0.6230 0.4851 -0.5774 -0.2417 0.4851 -0.5774 0.7439D = 1.0000 0 0 0 0.0000 00 0 1.0000贊同四、所得圖形kmax(y(k)x(k) = y(k)/max(x(k)x

14、(k+1) = Ay(k)01(1, 1, 1)(10, 8, 1)110(1, 0.8, 0.1)(7.2, 5.4, -0.8)27.2(1, 0.75, -0.111111)(6.5, 4.75, -1.222222)36.57(1, 0.730769, -0.203704)(6.230766, 4.499997, -1.407408)46.230766(1, 0.722222, -0.225880)(6.111108, 4.388886, -1.1451767)56.111108(1, 0.718182, -0.237561)(6.054548, 4.336336, -1.475122

15、)66.054548(1, 0.716216, -0.243639)(6.027024, 4.310808, -1.487278)76.027024(1, 0.715247, -0.246768)(6.013458, 4.298211, -1.483536)86.013458(1, 0.714765, -0.248366)(6.00671, 4.291945, -1.496732)96.00671(1, 0.714525, -0.249177)(6.00335, 4.28825, -1.496354)106.00335(1, 0.714405, -0.249586)(6.00167, 4.28

16、7265, -1.499172)116.00167(1, 0.714345, -0.239792)(6.00083, 4.286485, -1.499584)126.00083(1, 0.714315, -0.249896)五、經(jīng)驗總結：MATLAB提供了eig來計算矩陣的特征值、特征向量信息。如果再結合使用MATLAB的排序函數(shù)等資源，可以綜合利用來解決問題。實驗四 數(shù)值積分一、算法設計思想：辛普生公式,亦即為拋物線法求其近似值,其實質(zhì)是將曲線弧換為二次拋物線弧,再對拋物線積分。其結果為:由于“曲線”y=f(x)在區(qū)間a,b上的平均“高度(”值)為:(如圖)代入(1)式得:(2)其中:yi(

17、i=0,1,2,32n)為將區(qū)間a,b進行2n等份,得2n+1個分點,依次每個分點對應的函數(shù)值。龍貝格求積的基本思路與計算步驟：用復化梯形公式，取區(qū)間長度為，復化梯形公式的值與原積分值之間存在關系： 該式稱為歐拉-麥克勞林求和公式。上式還表明，用近似的截斷誤差為：龍貝格求積算法計算步驟如下：(1) 計算：(2) 對分區(qū)間并計算和：，其中為新分點的函數(shù)值之和。(3) 計算與：，(4) 利用外推公式：，()，()，直至求出。(5) 判斷是否真，其中為給定的精度。若真，則，否則重復(2)、(3)、(4)的計算。排成三角數(shù)表數(shù)表沿豎向和斜向都是收斂于，即當趨于無窮，與均收斂于。二、程序清單：第一題：用

18、梯形公式和辛浦生公式計算積分，并估計誤差。（計算取5位小數(shù)）int main() int a=1 ;int b=2 ;int n=100 ; double h = 0.01 ;double xk = 0 ;double xk1 = 0 ;double xk12 = 0 ;double i = 0 ;for ( int k =0; k<=n-1; k+ )xk = a + k*h ;xk12 = xk + h/2 ;xk1 = xk + h ;i = i + f(xk) + 4*f(xk12) + f( xk1 ) ; i = i * h/6 ; cout << "

19、the value is " << i << endl; cout << " Done ! " << endl; return 0 ;第二題：用龍貝格法求積分，要求誤差不超過。（計算取6位小數(shù)）functionR,quad,err,h=romber(f,a,b,n,tol)M=1;h=b-a;err=1;J=0;R=zeros(4,4);R(1,1)=h*(feval('f',a)+feval('f',b)/2;while(err>tol)&(J<n)|(J<4) J=J+1; h=h/2; s=0; for p=1:M x=a+h*(2*p-1); s=s+feval('f',x); end R(J+1,1)=R(J,1)/2+h*s; M=2*M; for K=1:J R(J+1,K+1)=R(J+1,K)+(R(J+