計(jì)算方法實(shí)驗(yàn)報(bào)告習(xí)題浙大

1、專業(yè)：電氣工程及其自動(dòng)化姓名： 李X 計(jì)算方法實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱： 實(shí)驗(yàn)1 從函數(shù)表出發(fā)進(jìn)行插值 1 引言某個(gè)實(shí)際問題中，函數(shù)f(x)在區(qū)間a,b上存在且連續(xù)，但難以找到其表達(dá)式，只能通過實(shí)驗(yàn)和觀測(cè)得到有限點(diǎn)上的函數(shù)表。有些情況雖然可以寫出表達(dá)式，但結(jié)構(gòu)復(fù)雜，使用不方便。所以希望構(gòu)造簡(jiǎn)單函數(shù)P(x)作為f(x)的近似值。插值法是解決此類問題的一種方法。設(shè)函數(shù)y=在插值區(qū)間a,b上連續(xù)，且在n+1個(gè)不同的插值節(jié)點(diǎn)ax0,x1,xnb上分別取值y0,y1,yn。目的是要在一個(gè)性質(zhì)優(yōu)良、便于計(jì)算的插值函數(shù)類中，求一簡(jiǎn)單函數(shù)P(x)，滿足插值條件P(xi)=yi(i=0,1,n)，而在其他點(diǎn)xxi上，作

2、為f(x)近似值。求插值函數(shù)P(x)的方法稱為插值法1。2 實(shí)驗(yàn)?zāi)康暮鸵筮\(yùn)用Matlab編寫m文件，定義三種插值函數(shù)，要求一次性輸入整張函數(shù)表，并利用計(jì)算機(jī)選擇在插值計(jì)算中所需的節(jié)點(diǎn)。分別通過分段線性插值、分段二次插值和全區(qū)間上拉格朗日插值計(jì)算f(0.15)，f(0.31)，f(0.47)的近似值。3 算法原理與流程圖（1）原理1. 線性插值當(dāng)給定了n+1個(gè)點(diǎn)x0<x1<<xn上的函數(shù)值y0,y1,yn后，若要計(jì)算xxi處函數(shù)值f(x)的近似值，可先選取兩個(gè)節(jié)點(diǎn)xi-1與xi使xxi-1,xi，然后在小區(qū)間xi-1,xi上作線性插值，即得這種分段低次插值叫分段線性插值。2.

3、 分段二次插值當(dāng)給定了n+1個(gè)點(diǎn)x0<x1<<xn上的函數(shù)值y0,y1,yn后，若要計(jì)算xxi處函數(shù)值f(x)的近似值，可先選取距離x最近的三個(gè)節(jié)點(diǎn)xi-1，xi與xi+1，然后進(jìn)行二次插值，即得這種分段低次插值叫分段二次插值。3. 全區(qū)間上拉格朗日插對(duì)節(jié)點(diǎn)xi(i=0,1,n)中任一點(diǎn)xk(0kn)，作一n次多項(xiàng)式lk(x)，使它在該點(diǎn)上的取值為1，在其余點(diǎn)xi(i=0,1,k-1,k+1,n)上取值為零。對(duì)應(yīng)于每一節(jié)點(diǎn)xk(k=0,1,n)，都能寫出一個(gè)滿足此條件的多項(xiàng)式，這樣寫出n+1個(gè)多項(xiàng)式l0(x),l1(x),ln(x)。拉格朗日n次插值多項(xiàng)式（對(duì)于全區(qū)間上的插值

4、，n取函數(shù)表的長(zhǎng)度）（2） 流程圖4 程序代碼及注釋1. 分段線性插值法%分段線性插值function y=fdxx(x0,y0,x) %定義函數(shù)p=length(y0);n=length(x0);m=length(x); %計(jì)算函數(shù)表和x的長(zhǎng)度if p=n error('數(shù)據(jù)輸入有誤，請(qǐng)重新輸入');%若函數(shù)表的x與y長(zhǎng)度不一致則輸入有誤else fprintf('分段線性插值nn'); for t=1:m %利用循環(huán)計(jì)算每個(gè)x的插值結(jié)果 z=x(t); if z<x0(1)|z>x0(n) fprintf('x(%d)超出范圍;n'

5、;,t); break; %若x不在函數(shù)表范圍內(nèi)，則插值結(jié)果將不準(zhǔn)確 end for i=1:n-1 if z<x0(i+1) break; %選取合適的兩點(diǎn)使x(i)<x<x(i+1) end end %注：若x不在函數(shù)表范圍內(nèi)，則i=n-1 y(t)=y0(i)*(z-x0(i+1)/(x0(i)-x0(i+1)+y0(i+1)*(z-x0(i)/(x0(i+1)-x0(i); %按照分段線性插值公式求解y fprintf('y(%d)=%fnx1=%.3f y1=%.3f,x2=%.3f y2=%.3fnn',t,y(t),x0(i),y0(i),x0(

6、i+1),y0(i+1); %輸出插值結(jié)果和所需的節(jié)點(diǎn) endend注：若要在Matlab中直接調(diào)用插值函數(shù)，需要事先將習(xí)題1函數(shù)表輸入2. 分段二次插值法%分段二次插值function y=fd2(x0,y0,x) %定義函數(shù)分段二次插值p=length(y0);n=length(x0);m=length(x); %計(jì)算函數(shù)表和x的長(zhǎng)度if p=n error('數(shù)據(jù)輸入有誤，請(qǐng)重新輸入');%若函數(shù)表的x與y長(zhǎng)度不一致則輸入有誤else fprintf('分段二次差值nn'); for t=1:m %運(yùn)用循環(huán)求解所有點(diǎn)的插值 z=x(t); if z<

7、x0(1)|z>x0(n) fprinf('x(%d)超出范圍;n',t); break; %如果x不在函數(shù)表范圍內(nèi)無法插值 end i=n-1; %若下列循環(huán)i不變，則i=n-1， for j=1:n-2 if z<(x0(j+1)+x0(j)/2 i=j; end end %選取i使得x(i-1),x(i),x(i+1)是距x最近的三個(gè)點(diǎn) s=0.0; for k=i-1:i+1 p=1; for l=i-1:i+1 if l=k p=p*(z-x0(l)/(x0(k)-x0(2); end end s=s+y0(k)*p; end %根據(jù)分段二次插值公式求y

8、y(t)=s; fprintf('y(%d)=%fnx1=%.3f y1=%.3fnx2=%.3f y2=%.3fnx3=%.3f y3=%.3fnn',t,y(t),x0(i-1),y0(i-1),x0(i),y0(i),x0(i+1),y0(i+1); %輸出結(jié)果和所需節(jié)點(diǎn) endend3.拉格朗日全區(qū)間上插值%拉格朗日插值function y=lagrange(x0,y0,x) %定義函數(shù)p=length(y0);n=length(x0);m=length(x); %計(jì)算函數(shù)表和x的長(zhǎng)度if p=n error('數(shù)據(jù)輸入有誤，請(qǐng)重新輸入');%若函數(shù)表的

9、x與y長(zhǎng)度不一致則輸入有誤else fprintf('拉格朗日插值n'); for t=1:m %利用循環(huán)計(jì)算每個(gè)x的插值 s=0.0; z=x(t); for k=1:n p=1; for i=1:n if i=k p=p*(z-x0(j)/(x0(k)-x0(i); end end s=s+y0(k)*p; end %根據(jù)拉格朗日插值公式求解y y(t)=s; fprintf('y(%d)=%fn',t,y(t); %輸出插值結(jié)果 endend5 算例分析x0.00.10.1950.30.4010.5f(x)0.398940.396950.391420.38

10、1380.368120.35206從函數(shù)表出發(fā)計(jì)算f(0.15)，f(0.31)，f(0.47)的近似值>> x0=0.0 0.1 0.195 0.3 0.401 0.5;>> y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206;>> x=0.15 0.31 0.47;%首先在Matlab中輸入函數(shù)表和待插值變量x%確保Current Directory 中存放分別含有上述代碼的三個(gè)m文件（1） 分段線性插值>> y1=fdxx(x0,y0,x) %調(diào)用分段線性插值函數(shù)得到下列結(jié)果分段線性插值y(

11、1)=0.394039x1=0.100 y1=0.397,x2=0.195 y2=0.391y(2)=0.380067x1=0.300 y1=0.381,x2=0.401 y2=0.368y(3)=0.356927x1=0.401 y1=0.368,x2=0.500 y2=0.352y1 = 0.39404 0.38107 0.35693（2） 分段二次插值>> y=fd2(x0,y0,x) %調(diào)用分段二次插值函數(shù)得到下列結(jié)果分段二次差值y(1)=0.394554x1=0.195 y1=0.391x2=0.300 y2=0.381x3=0.401 y3=0.368y(2)=0.38

12、0225x1=0.195 y1=0.391x2=0.300 y2=0.381x3=0.401 y3=0.368y(3)=0.357247x1=0.300 y1=0.381x2=0.401 y2=0.368x3=0.500 y3=0.352y = 0.39455 0.38122 0.35725（3） 全區(qū)間上拉格朗日插值>> y=lagrange(x0,y0,x) %調(diào)用拉格朗日性插值函數(shù)得到下列結(jié)果拉格朗日插值y(1)=0.394473y(2)=0.380219y(3)=0.357222y = 0.39447 0.38122 0.35722（4） 測(cè)試示例>> x=0.

13、15 0.31 0.51;>> y=fdxx(x0,y0,x)分段線性插值y(1)=0.394039x1=0.100 y1=0.397,x2=0.195 y2=0.391y(2)=0.380067x1=0.300 y1=0.381,x2=0.401 y2=0.368x(3)超出范圍;>> y1= 0.39894 0.39695 0.39142 0.38138 0.36812 ;>> y=fdxx(x0,y1,x)? Error using => fdxx數(shù)據(jù)輸入有誤，請(qǐng)重新輸入6 討論與結(jié)論1. 對(duì)程序運(yùn)行的時(shí)間加以比較使用tic,toc函數(shù)計(jì)算下列四

14、種方法計(jì)算上述問題所運(yùn)行的時(shí)間函數(shù)interp1(x0,y0,x)fdxx(x0,y0,x)fd2(x0,y0,x)lagrange(x0,y0,x)運(yùn)行時(shí)間(s)0.296990.0629990.0764340.0348170.255580.0723150.0718960.0382890.242970.0699530.0850470.048222從三次實(shí)驗(yàn)結(jié)果可知，interp1函數(shù)所用時(shí)間比編寫的fdxx函數(shù)長(zhǎng)（interp1函數(shù)為matlab求分段線性插值的函數(shù)）；而fdxx和另外兩個(gè)函數(shù)時(shí)間比較接近2. 程序優(yōu)化由分段線性插值和分段二次插值的原理，x取值在函數(shù)表范圍內(nèi)時(shí)，插值結(jié)果有意義，而當(dāng)x取值在函數(shù)表范圍以外，利用分段線性插值公式仍可以進(jìn)行運(yùn)算并得到一個(gè)值，但其結(jié)果不準(zhǔn)確；分段二次插值則無法找到三個(gè)合適的點(diǎn)以求插值，不予以輸出結(jié)果；若輸入的函數(shù)表x與y的長(zhǎng)度不相等，則無法插值。所以加入判斷p=length(y0);n=length(x0);if p=n error('數(shù)據(jù)輸入有誤，請(qǐng)重新輸入');if z<x0(1)|z>x0(n) fprintf('x(%d)超出范圍;n',t); break;end從而避免了因輸入錯(cuò)誤導(dǎo)致