數(shù)學建模競賽常用方法之數(shù)據(jù)處理

數(shù)學建模競賽常用方法之數(shù)據(jù)處理方法南京師范大學泰州學院數(shù)科院王明剛主要內(nèi)容一、

數(shù)據(jù)的導入與導出三、

數(shù)據(jù)插值二、

數(shù)據(jù)預處理四、

數(shù)據(jù)擬合2024/8/3從TXT文件中讀取數(shù)據(jù)MATLAB中讀取文本文件的常用函數(shù)2024/8/3一、利用數(shù)據(jù)導入向導導入TXT文件2024/8/3數(shù)據(jù)文件略去……【例2.1-1】利用數(shù)據(jù)導入向導讀取文件examp02_01.txt至examp02_11.txt中的數(shù)據(jù)2024/8/3二、調(diào)用高級函數(shù)讀取數(shù)據(jù)1.調(diào)用importdata函數(shù)讀取數(shù)據(jù)調(diào)用格式：importdata(filename)A=importdata(filename)A=importdata(filename,delimiter)A=importdata(filename,delimiter,headerline)[AD]=importdata(…)[ADH]=importdata(…)[…]=importdata(‘-pastespecial’,…)

2024/8/3>>importdata('examp02_04.txt')>>x=importdata('examp02_07.txt')>>x=importdata('examp02_03.txt',';')>>x=importdata('examp02_08.txt','',2)>>[x,s,h]=importdata('examp02_07.txt')>>FileContent=importdata('examp02_10.txt')【例2.1-2】調(diào)用importdata函數(shù)讀取文件examp02_01.txt至examp02_11.txt中的數(shù)據(jù)2024/8/32.調(diào)用load函數(shù)讀取數(shù)據(jù)調(diào)用格式：S=load(filename)S=load(filename,variables)S=load(filename,'-mat',variables)S=load(filename,'-ascii')load(...)load...2024/8/3>>loadexamp02_01.txt>>load-asciiexamp02_01.txt>>x1=load('examp02_02.txt')>>x1=load('examp02_02.txt','-ascii');>>loadexamp02_03.txt>>loadexamp02_04.txt……【例2.1-3】調(diào)用load函數(shù)讀取文件examp02_01.txt至examp02_12.txt中的數(shù)據(jù)2024/8/33.調(diào)用dlmread函數(shù)讀取數(shù)據(jù)調(diào)用格式：M=dlmread(filename)M=dlmread(filename,delimiter)M=dlmread(filename,delimiter,R,C)M=dlmread(filename,delimiter,range)2024/8/3>>x=dlmread('examp02_03.txt')>>x=dlmread('examp02_03.txt',',',2,3)>>x=dlmread('examp02_03.txt',',',[1,2,2,5])>>x=dlmread('examp02_05.txt')>>x=dlmread('examp02_06.txt')>>x=dlmread('examp02_09.txt')……【例2.1-4】調(diào)用dlmread函數(shù)讀取文件examp02_01.txt至examp02_11.txt中的數(shù)據(jù)2024/8/34.調(diào)用textread函數(shù)讀取數(shù)據(jù)調(diào)用格式：[A,B,C,…]=textread('filename','format')[A,B,C,…]=textread('filename','format',N)[…]=textread(…,'param','value',…)textread函數(shù)支持的format字符串textread函數(shù)支持的參數(shù)名與參數(shù)值列表2024/8/3>>x1=textread('examp02_01.txt');>>x2=textread('examp02_02.txt');>>x3=textread('examp02_03.txt','','delimiter',',');>>[c1,c2,c3,c4,c5]=textread('examp02_04.txt','%f%f%f

%f%f','delimiter',',;*');>>x5=textread('examp02_05.txt','','emptyvalue',-1)>>x8=textread('examp02_08.txt','','headerlines',7)>>x9=textread('examp02_09.txt','','delimiter',',','whitespace','+i')>>[c1,c2,c3,c4,c5,c6,c7,c8]=textread('examp02_09.txt',...'%f%f%f

%f

%f

%f

%f%f','delimiter',',','whitespace','+i');……【例2.1-5】調(diào)用textread函數(shù)讀取文件examp02_01.txt至examp02_11.txt中的數(shù)據(jù)2024/8/3三、調(diào)用低級函數(shù)讀取數(shù)據(jù)1.調(diào)用fopen函數(shù)打開文件調(diào)用格式：[fid,message]=fopen(filename,permission)[filename,permission]=fopen(fid)2024/8/32.調(diào)用fclose函數(shù)關閉文件調(diào)用格式：status=fclose(fid)status=fclose('all')2024/8/33.調(diào)用fseek、ftell、frewind和feof函數(shù)控制讀寫位置調(diào)用格式：status=fseek(fid,offset,origin)%設定文件指針位置position=ftell(fid)%獲取文件指針位置frewind(fid)%移動當前文件指針到文件的開頭eofstat=feof(fid)%判斷是否到達文件末尾2024/8/34.調(diào)用fgets、fgetl函數(shù)讀取文件的下一行調(diào)用格式：tline=fgets(fid)tline=fgets(fid,nchar)tline=fgetl(fid)2024/8/35.調(diào)用textscan函數(shù)讀取數(shù)據(jù)調(diào)用格式：C=textscan(fid,'format')C=textscan(fid,'format',N)C=textscan(fid,'format',param,value,…)C=textscan(fid,'format',N,param,value,…)C=textscan(str,…)[C,position]=textscan(…)2024/8/3>>fid=fopen('examp02_08.txt','r');%以只讀方式打開文件examp02_08.txt>>fgets(fid);%讀取文件的第1行>>fgets(fid);%讀取文件的第2行>>A=textscan(fid,'%f%f%f

%f

%f

%f','CollectOutput',1)A=[3x6double]>>fgets(fid);%讀取文件的第6行>>fgets(fid);%讀取文件的第7行>>B=textscan(fid,'%f%f%f','CollectOutput',1)B=[2x3double]>>fclose(fid);%關閉文件……【例2.1-6】調(diào)用textscan函數(shù)讀取文件examp02_01.txt至examp02_13.txt中的數(shù)據(jù)2024/8/3把數(shù)據(jù)寫入TXT文件MATLAB中寫文本文件的常用函數(shù)2024/8/3一、調(diào)用dlmwrite函數(shù)寫入數(shù)據(jù)調(diào)用格式：dlmwrite(filename,M)dlmwrite(filename,M,'D')dlmwrite(filename,M,'D',R,C)dlmwrite(filename,M,'attrib1',value1,'attrib2',value2,…)dlmwrite(filename,M,'-append')dlmwrite(filename,M,'-append',attribute-valuelist)dlmwrite函數(shù)支持的參數(shù)名與參數(shù)值列表2024/8/3%定義復數(shù)矩陣>>x=[1.455390+1.360686i8.692922+5.797046i5.498602+1.449548i8.530311+6.220551i3.509524+5.132495i4.018080+0.759667i2.399162+1.233189i1.839078+2.399525i4.172671+0.496544i9.027161+9.447872i4.908641+4.892526i3.377194+9.000538i];%將復數(shù)矩陣x寫入文件examp02_09.txt，用逗號(',')作分隔符，用'\r\n'作換行符>>dlmwrite('examp02_09.txt',x,'delimiter',',','newline','pc')【例2.2-1】用逗號作為分隔符，調(diào)用dlmwrite函數(shù)將如下復數(shù)矩陣寫入文件examp02_09.txt2024/8/3二、調(diào)用fprintf函數(shù)寫入數(shù)據(jù)調(diào)用格式：count=fprintf(fid,format,A,…)【例】

%在屏幕上顯示一句話>>y=fprintf(1,'祝福我們偉大的新中國%d周歲生日快樂?。?！',60)祝福我們偉大的新中國60周歲生日快樂?。?！y=382024/8/3%產(chǎn)生一個8行5列的隨機矩陣，其元素服從[0,10]上的均勻分布>>x=10*rand(8,5);%以寫入方式打開文件，返回文件標識符>>fid=fopen('examp02_01.txt','wt');%把矩陣x以指定格式寫入文件examp02_01.txt>>fprintf(fid,'%-f%-f%-f

%-f

%-f

%-f

%-f%-f\n',x);>>fclose(fid);%關閉文件……【例2.2-2】用fprintf函數(shù)將數(shù)據(jù)寫入文件examp02_01.txt至examp02_11.txt的代碼調(diào)用fprintf函數(shù)寫入數(shù)據(jù)或在屏幕上顯示數(shù)據(jù)時，format參數(shù)指定的格式循環(huán)作用在矩陣的列上，原始矩陣的列在文件中或屏幕上就變成了行。2024/8/3從Excel文件中讀取數(shù)據(jù)2024/8/3一、利用數(shù)據(jù)導入向導導入Excel文件2024/8/3二、調(diào)用xlsread函數(shù)讀取數(shù)據(jù)調(diào)用格式：num=xlsread(filename)num=xlsread(filename,-1)num=xlsread(filename,sheet)num=xlsread(filename,range)num=xlsread(filename,sheet,range)num=xlsread(filename,sheet,range,'basic')num=xlsread(filename,…,functionhandle)[num,txt]=xlsread(filename,…)[num,txt,raw]=xlsread(filename,…)[num,txt,raw,X]=xlsread(filename,…,functionhandle)xlsreadfilenamesheetrangebasic2024/8/3%讀取文件examp02_14.xls第1個工作表中單元格A2:H4中的數(shù)據(jù)>>num=xlsread('examp02_14.xls','A2:H4')%返回讀取的數(shù)據(jù)矩陣numnum=1601016010101NaN063632601016010102NaN073733601016010103NaN000【例2.3-1】用xlsread函數(shù)的這種調(diào)用讀取文件examp02_14.xls第1個工作表中區(qū)域A2:H4的數(shù)據(jù)2024/8/3把數(shù)據(jù)寫入Excel文件2024/8/3一、調(diào)用xlswrite函數(shù)寫數(shù)據(jù)到Excel調(diào)用格式：xlswrite(filename,M)xlswrite(filename,M,sheet)xlswrite(filename,M,range)xlswrite(filename,M,sheet,range)status=xlswrite(filename,…)[status,message]=xlswrite(filename,…)xlswritefilenameMsheetrange2024/8/3%生成一個10行10列的隨機矩陣，其元素服從[0,1]上的均勻分布>>x=rand(10);%把矩陣x寫入文件examp02_15.xls的第2個工作表中的單元格區(qū)域D6:M15，并返回操作信息>>[s,t]=xlswrite('examp02_15.xls',x,2,'D6:M15')s=1t=message:''identifier:''【例2.4-1】生成一個10×10的隨機數(shù)矩陣，將它寫入Excel文件examp02_15.xls的第2個工作表的指定區(qū)域2024/8/3>>x={1,60101,6010101,'陳亮',63,'';2,60101,6010102,'李旭',73,'';3,60101,...6010103,'劉鵬飛',0,'缺考'}%定義一個元胞數(shù)組x=[1][60101][6010101]'陳亮'[63]''[2][60101][6010102]'李旭'[73]''[3][60101][6010103]'劉鵬飛'[0]'缺考'%把元胞數(shù)組x寫入文件examp02_15.xls的指定工作表（xiezhh）中的單元格區(qū)域A3:F5>>xlswrite('examp02_15.xls',x,'xiezhh','A3:F5')【例2.4-2】定義一個元胞數(shù)組，將它寫入Excel文件examp02_15.xls的自命名工作表的指定區(qū)域2024/8/3%讀取文件examp02_14.xls中的數(shù)值型數(shù)據(jù)>>num=xlsread('examp02_14.xls');%把讀取的數(shù)據(jù)num寫入文件examp02_15.xls中的名稱為'成績'的工作表的默認區(qū)域>>xlswrite('examp02_15.xls',num,'成績')【例2.4-3】讀取文件examp02_14.xls中的數(shù)值型數(shù)據(jù)，并將讀取的數(shù)據(jù)寫入文件examp02_15.xls的自命名工作表的默認區(qū)域2024/8/3主要內(nèi)容數(shù)據(jù)的平滑處理數(shù)據(jù)的標準化變換數(shù)據(jù)的極差歸一化變換數(shù)據(jù)的預處理2024/8/3數(shù)據(jù)的平滑處理2024/8/3一、smooth函數(shù)調(diào)用格式：yy=smooth(y)yy=smooth(y,span)yy=smooth(y,method)yy=smooth(y,span,method)yy=smooth(y,'sgolay',degree)yy=smooth(y,span,'sgolay',degree)yy=smooth(x,y,…)2024/8/3%產(chǎn)生一個從0到2*pi的向量，長度為500>>t=linspace(0,2*pi,500)';>>y=100*sin(t);%產(chǎn)生正弦波信號%產(chǎn)生500行1列的服從N(0,152)分布的隨機數(shù)，作為噪聲信號>>noise=normrnd(0,15,500,1);>>y=y+noise;%將正弦波信號加入噪聲信號>>figure;%新建一個圖形窗口>>plot(t,y);%繪制加噪波形圖>>xlabel('t');%為X軸加標簽>>ylabel('y=sin(t)+噪聲');%為Y軸加標簽【例3.1-1】產(chǎn)生一列正弦波信號，加入噪聲信號，然后調(diào)用smooth函數(shù)對加入噪聲的正弦波進行濾波（平滑處理）2024/8/3移動平均法：>>yy1=smooth(y,30);%利用移動平均法對y進行平滑處理>>figure;%新建一個圖形窗口>>plot(t,y,'k:');%繪制加噪波形圖>>holdon;>>plot(t,yy1,'k','linewidth',3);%繪制平滑后波形圖>>xlabel('t');%為X軸加標簽>>ylabel('moving');%為Y軸加標簽>>legend('加噪波形','平滑后波形');2024/8/3lowess方法：%利用lowess方法對y進行平滑處理>>yy2=smooth(y,30,'lowess');>>figure;%新建一個圖形窗口>>plot(t,y,'k:');%繪制加噪波形圖>>holdon;>>plot(t,yy2,'k','linewidth',3);%繪制平滑后波形圖>>xlabel('t');%為X軸加標簽>>ylabel('lowess');%為Y軸加標簽>>legend('加噪波形','平滑后波形');2024/8/3rlowess方法：%利用rlowess方法對y進行平滑處理>>yy3=smooth(y,30,'rlowess');>>figure;%新建一個圖形窗口>>plot(t,y,'k:');%繪制加噪波形圖>>holdon;>>plot(t,yy3,'k','linewidth',3);%繪制平滑后波形圖>>xlabel('t');%為X軸加標簽>>ylabel('rlowess');%為Y軸加標簽>>legend('加噪波形','平滑后波形');2024/8/3loess方法：%利用loess方法對y進行平滑處理>>yy4=smooth(y,30,'loess');>>figure;%新建一個圖形窗口>>plot(t,y,'k:');%繪制加噪波形圖>>holdon;>>plot(t,yy4,'k','linewidth',3);%繪制平滑后波形圖>>xlabel('t');%為X軸加標簽>>ylabel('loess');%為Y軸加標簽>>legend('加噪波形','平滑后波形');2024/8/3sgolay方法：%利用sgolay方法對y進行平滑處理>>yy5=smooth(y,30,'sgolay',3);>>figure;%新建一個圖形窗口>>plot(t,y,'k:');%繪制加噪波形圖>>holdon;>>plot(t,yy5,'k','linewidth',3);%繪制平滑后波形圖>>xlabel('t');%為X軸加標簽>>ylabel('sgolay');%為Y軸加標簽>>legend('加噪波形','平滑后波形');2024/8/32024/8/3二、smoothts函數(shù)調(diào)用格式：output=smoothts(input)output=smoothts(input,‘b’,wsize)%盒子法output=smoothts(input,‘g’,wsize,stdev)%高斯窗方法output=smoothts(input,‘e’,n)%指數(shù)法2024/8/3【例3.1-2】現(xiàn)有上海股市日開盤價、最高價、最低價、收盤價、收益率等數(shù)據(jù)，時間跨度為2005年1月4日至2007年4月3日，共510組數(shù)據(jù)。完整數(shù)據(jù)保存在文件examp03_02.xls中，其中部分數(shù)據(jù)如下圖所示。試調(diào)用smoothts函數(shù)對日收盤價數(shù)據(jù)進行平滑處理2024/8/3繪制日收盤價曲線圖：%從文件examp03_02.xls中讀取數(shù)據(jù)>>x=xlsread('examp03_02.xls');>>price=x(:,4)';%提取矩陣x的第4列數(shù)據(jù)，即收盤價數(shù)據(jù)>>figure;%新建一個圖形窗口%繪制日收盤價曲線圖，黑色實線，線寬為2>>plot(price,'k','LineWidth',2);%為X軸和Y軸加標簽>>xlabel('觀測序號');>>ylabel('上海股市日收盤價');2024/8/3盒子法：%用盒子法平滑數(shù)據(jù)，窗寬為30>>output1=smoothts(price,'b',30);%用盒子法平滑數(shù)據(jù)，窗寬為100>>output2=smoothts(price,'b',100);>>figure;%新建一個圖形窗口>>plot(price,'.');%繪制日收盤價散點圖>>holdon%繪制平滑后曲線圖，黑色實線，線寬為2>>plot(output1,'k','LineWidth',2);%繪制平滑后曲線圖，黑色點劃線，線寬為2>>plot(output2,'k-.','LineWidth',2);>>xlabel('觀測序號');ylabel('Boxmethod');%為X軸和Y軸加標簽%為圖形加標注框>>legend('原始散點','平滑曲線(窗寬30)','平滑曲線(窗寬100)','location','northwest');2024/8/3用高斯窗方法平滑數(shù)據(jù)：>>output3=smoothts(price,'g',30);%窗寬為30，標準差為默認值0.65>>output4=smoothts(price,'g',100,100);%窗寬為100，標準差為100>>figure;%新建一個圖形窗口>>plot(price,'.');%繪制日收盤價散點圖>>holdon%繪制平滑后曲線圖，黑色實線，線寬為2>>plot(output3,'k','LineWidth',2);%繪制平滑后曲線圖，黑色點劃線，線寬為2>>plot(output4,'k-.','LineWidth',2);%為X軸和Y軸加標簽>>xlabel('觀測序號');ylabel('Gaussianwindowmethod');>>legend('原始散點','平滑曲線(窗寬30，標準差0.65)',...'平滑曲線(窗寬100，標準差100)','location','northwest');2024/8/3用指數(shù)法平滑數(shù)據(jù)：>>output5=smoothts(price,'e',30);%用指數(shù)法平滑數(shù)據(jù)，窗寬為30>>output6=smoothts(price,'e',100);%用指數(shù)法平滑數(shù)據(jù)，窗寬為100>>figure;%新建一個圖形窗口>>plot(price,'.');%繪制日收盤價散點圖>>holdon%繪制平滑后曲線圖，黑色實線，線寬為2>>plot(output5,'k','LineWidth',2);%繪制平滑后曲線圖，黑色點劃線，線寬為2>>plot(output6,'k-.','LineWidth',2);%為X軸和Y軸加標簽>>xlabel('觀測序號');ylabel('Exponentialmethod');>>legend('原始散點','平滑曲線(窗寬30)','平滑曲線(窗寬100)','location','northwest');2024/8/32024/8/3三、medfilt1函數(shù)（一維中值濾波）調(diào)用格式：y=medfilt1(x,n)y=medfilt1(x,n,blksz)y=medfilt1(x,n,blksz,dim)2024/8/3%產(chǎn)生一個從0到2*pi的向量，長度為500>>t=linspace(0,2*pi,500)';>>y=100*sin(t);%產(chǎn)生正弦波信號%產(chǎn)生500行1列的服從N(0,152)分布的隨機數(shù)，作為噪聲信號>>noise=normrnd(0,15,500,1);>>y=y+noise;%將正弦波信號加入噪聲信號>>figure;%新建一個圖形窗口>>plot(t,y);%繪制加噪波形圖>>xlabel('t');%為X軸加標簽>>ylabel('y=sin(t)+噪聲');%為Y軸加標簽【例3.1-3】產(chǎn)生一列正弦波信號，加入噪聲信號，然后調(diào)用medfilt1函數(shù)對加入噪聲的正弦波進行濾波（平滑處理）2024/8/3中值濾波：%調(diào)用medfilt1對加噪正弦波信號y進行中值濾波，并繪制波形圖>>yy=medfilt1(y,30);%指定窗寬為30，對y進行中值濾波>>figure;%新建一個圖形窗口>>plot(t,y,'k:');%繪制加噪波形圖>>holdon%繪制平滑后曲線圖，黑色實線，線寬為3>>plot(t,yy,'k','LineWidth',3);>>xlabel('t');%為X軸加標簽>>ylabel('中值濾波');%為Y軸加標簽>>legend('加噪波形','平滑后波形');2024/8/32024/8/3數(shù)據(jù)的標準化變換2024/8/3一、標準化變換公式2024/8/3二、標準化變換的MATLAB實現(xiàn)調(diào)用格式：Z=zscore(X)[Z,mu,sigma]=zscore(X)[…]=zscore(X,1)[…]=zscore(X,flag,dim)1.zscore函數(shù)2024/8/3%調(diào)用rand函數(shù)產(chǎn)生一個10行，4列的隨機矩陣，每列服從不同的均勻分布>>x=[rand(10,1),5*rand(10,1),10*rand(10,1),500*rand(10,1)]%調(diào)用zscore函數(shù)對x進行標準化變換（按列標準化），%返回變換后矩陣xz，以及矩陣x各列的均值構成的向量mu，各列的標準差構成的向量sigma>>[xz,mu,sigma]=zscore(x)>>std(xz)%求標準化后矩陣xz的各列的標準差【例3.2-1】調(diào)用rand函數(shù)產(chǎn)生一個隨機矩陣，然后調(diào)用zscore函數(shù)將其按列標準化2024/8/3數(shù)據(jù)的極差歸一化變換2024/8/3一、極差歸一化變換公式2024/8/3二、極差歸一化變換的MATLAB實現(xiàn)調(diào)用格式：R=rscore(X)[R,xmin,xrange]=rscore(X)[…]=rscore(X,dim)1.自編的rscore函數(shù)2024/8/3%調(diào)用rand函數(shù)產(chǎn)生一個10行，4列的隨機矩陣，每列服從不同的均勻分布>>x=[rand(10,1),5*rand(10,1),10*rand(10,1),500*rand(10,1)]%調(diào)用rscore函數(shù)對x按列進行極差歸一化變換，%返回變換后矩陣R，以及矩陣x各列的最小值構成的向量xmin，各列的極差構成的向量xrange>>[R,xmin,xrange]=rscore(x)【例3.3-1】調(diào)用rand函數(shù)產(chǎn)生一個隨機矩陣，然后調(diào)用rscore函數(shù)對其按列進行極差歸一化變換2024/8/3調(diào)用格式：[Y,PS]=mapminmax(X,YMIN,YMAX)[Y,PS]=mapminmax(X,FP)Y=mapminmax('apply',X,PS)X=mapminmax('reverse',Y,PS)dy_dx=mapminmax('dy_dx',X,Y,PS)dx_dy=mapminmax('dx_dy',X,Y,PS)2.自帶的mapminmax函數(shù)2024/8/3%調(diào)用rand函數(shù)產(chǎn)生一個10行，4列的隨機矩陣，每列服從不同的均勻分布>>x=[rand(10,1),5*rand(10,1),10*rand(10,1),500*rand(10,1)]%調(diào)用rscore函數(shù)對x按列進行極差歸一化變換，%返回變換后矩陣R，以及矩陣x各列的最小值構成的向量xmin，各列的極差構成的向量xrange>>FP.ymin=0;>>FP.ymax=1;>>y=mapminmax(x',FP)'【例3.3-2】調(diào)用rand函數(shù)產(chǎn)生一個隨機矩陣，然后調(diào)用mapminmax函數(shù)對其按列進行極差歸一化變換數(shù)據(jù)插值基本定義：插值問題：給定一批數(shù)據(jù)點，需確定滿足特定要求的曲線或曲面，如果要求所求曲線（面）通過所有給定的數(shù)據(jù)點，這就是插值問題。數(shù)據(jù)擬合：如果不要求曲線（面）通過所有給定的數(shù)據(jù)點，而是要求它反映對象整體的變化趨勢，可得到更簡單實用的近似函數(shù)，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。一維插值的定義已知n+1個節(jié)點其中互不相同，不妨設求任一插值點處的插值

節(jié)點可視為由產(chǎn)生,表達式復雜,或無封閉形式,或未知.

構造一個(相對簡單的)函數(shù)通過全部節(jié)點,即再用計算插值，即

返回1.拉格朗日(Lagrange)插值兩點一次(線性)插值多項式(一次拉格朗日插值公式):例1根據(jù)下表給出的平方根值，用線性插值計算sqrt(5)x14916Sqrt(x)1234三點二次(拋物)插值公式:例2

根據(jù)例1給出的數(shù)據(jù)，用拋物線法計算x14916Sqrt(x)1234稱為拉格朗日插值基函數(shù)．

已知函數(shù)f(x)在n+1個點x0,x1,…,xn處的函數(shù)值為y0,y1,…,yn

．求一n次多項式函數(shù)Pn(x)，使其滿足：

Pn(xi)=yi,i=0,1,…,n.解決此問題的拉格朗日插值多項式公式如下其中Li(x)為n次多項式：一般情形采用拉格朗日多項式插值：選取不同插值節(jié)點n+1個，其中n為插值多項式的次數(shù)，當n分別取2,4,6,8,10時，繪出插值結果圖形.例1clcx=-5:0.01:5;y=1./(1+x.^2);plot(x,y,'linewidth',2)gtext('1/(1+x^2)')holdonpausen=3x0=linspace(-5,5,3);y0=1./(1+x0.^2);y1=lagr1(x0,y0,x)plot(x,y1,'r','linewidth',2)gtext('n=2')pausen=5x0=linspace(-5,5,5);y0=1./(1+x0.^2);y2=lagr1(x0,y0,x)plot(x,y2,'b','linewidth',2)gtext('n=4')pause拉格朗日多項式插值的這種振蕩現(xiàn)象叫Runge現(xiàn)象2.分段線性插值計算量與n無關;n越大，誤差越小.

xjxj-1xj+1x0xnxＯy例2用分段線性插值法求插值,并觀察插值誤差.1.在[-6,6]中平均選取5個點作插值4.在[-6,6]中平均選取41個點作插值2.在[-6,6]中平均選取11個點作插值3.在[-6,6]中平均選取21個點作插值clcx=-6:0.1:6;y=1./(1+x.^2);plot(x,y,'linewidth',2)gtext('1/(1+x^2)')holdonpausen=5x1=linspace(-6,6,5)y1=1./(1+x1.^2);plot(x1,y1,'linewidth',2)gtext('n=4')pausen=11x2=linspace(-6,6,11)y2=1./(1+x2.^2);plot(x2,y2,'linewidth',2)gtext('n=10')pause比分段線性插值更光滑

xyxi-1

xiab在數(shù)學上，光滑程度的定量描述是：函數(shù)(曲線)的k階導數(shù)存在且連續(xù)，則稱該曲線具有k階光滑性．光滑性的階次越高，則越光滑．是否存在較低次的分段多項式達到較高階光滑性的方法？三次樣條插值就是一個很好的例子．3.三次樣條插值

三次樣條插值g(x)為被插值函數(shù)．)()(limxgxSn=￥?例3用三次樣條插值選取11個基點計算插值clcx=-6:0.1:6;y=1./(1+x.^2);x1=linspace(-6,6,11);y1=interp1(x,y,x1,'spline')plot(x,y,'k',x1,y1,'r')用MATLAB作插值計算一維插值函數(shù)：yi=interp1(x，y，xi，'method')插值方法被插值點插值節(jié)點xi處的插值結果‘nearest’

最鄰近插值；‘linear’

線性插值；‘spline’

三次樣條插值；‘cubic’

立方插值；缺省時分段線性插值．注意：所有的插值方法都要求x是單調(diào)的，并且xi不能夠超過x的范圍．例4：從1點12點內(nèi)，每隔1小時測量一次溫度，測得的溫度的數(shù)值依次為：5，8，9，15，25，29，31，30，22，25，27，24．(1)試估計在3.1h,6.8h,7.5h,11.8h時的溫度值；(2)試估計每隔1/10小時的溫度值．clcsj=1:12;wd=[5,8,9,15,25,29,31,30,22,25,27,24];%sj1=[3.16.87.511.8];%t=interp1(sj,wd,sj1)%T=interp1(sj,wd,sj1,'spline')sj2=1:0.1:12;t=interp1(sj,wd,sj2,'spline')plot(sj,wd,'o',sj2,t,'*')例5已知飛機下輪廓線上數(shù)據(jù)如下，求x每改變0.1時的y值．

已知m

n個節(jié)點其中互不相同，不妨設構造一個二元函數(shù)通過全部已知節(jié)點,即再用計算插值，即;二維插值的定義

xyO第一種（網(wǎng)格節(jié)點）：第二種（散亂節(jié)點）：

yxO注意：最鄰近插值一般不連續(xù)．具有連續(xù)性的最簡單的插值是分片線性插值．1.最鄰近插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O二維或高維情形的最鄰近插值，與被插值點最鄰近的節(jié)點的函數(shù)值即為所求．網(wǎng)格節(jié)點插值有如下幾種形式：將四個插值點（矩形的四個頂點）處的函數(shù)值依次簡記為：

2.分片線性插值xy

(xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4插值函數(shù)為：第二片(上三角形區(qū)域)：(x,y)滿足插值函數(shù)為：注意：(x,y)當然應該是在插值節(jié)點所形成的矩形區(qū)域內(nèi)．顯然，分片線性插值函數(shù)是連續(xù)的；分兩片的函數(shù)表達式如下：第一片(下三角形區(qū)域)：(x,y)滿足返回雙線性插值是一片一片的空間二次曲面構成．雙線性插值函數(shù)的形式如下：其中有四個待定系數(shù)，利用該函數(shù)在矩形的四個頂點（插值節(jié)點）的函數(shù)值，得到四個代數(shù)方程，正好確定四個系數(shù)．3.雙線性插值x

y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O

要求x0,y0單調(diào)；x，y可取為矩陣，或x取行向量，y取為列向量，x,y的值分別不能超出x0,y0的范圍．z=interp2(x0,y0,z0,x,y,’method’)被插值點插值方法用MATLAB作網(wǎng)格節(jié)點數(shù)據(jù)的插值插值節(jié)點被插值點的函數(shù)值‘nearest’

最鄰近插值；‘linear’

雙線性插值；‘cubic’

雙三次插值；缺省時

雙線性插值.例：測得平板表面3×5網(wǎng)格點處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形．1.先在三維坐標畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲線圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個單位的地方進行插值.[x,y]=meshgrid(1:5,1:3);z=[8281808284;7963616581;8484828586];mesh(x,y,z)xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,z,xi',yi)mesh(xi,yi,zi)通過此例對最近鄰點插值、雙線性插值方法和雙三次插值方法的插值效果進行比較．x=1200:400:4000;y=1200:400:3600;z=xlsread('moutain.xls')mesh(x,y,z)figurexi=1200:50:4000;yi=1200:50:3600;zi1=interp2(x,y,z,xi',yi,'nearest')mesh(xi,yi,zi1)figurexi=1200:50:4000;yi=1200:50:3600;zi2=interp2(x,y,z,xi',yi)mesh(xi,yi,zi2)figurexi=1200:50:4000;yi=1200:50:3600;zi3=interp2(x,y,z,xi',yi,'cubic')mesh(xi,yi,zi3)figuresubplot(1,3,1)contour(xi,yi,zi1,1)

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點數(shù)據(jù)的插值計算

要求cx取行向量，cy取為列向量．被插值點插值方法插值節(jié)點被插值點的函數(shù)值‘nearest’最鄰近插值‘linear’

雙線性插值‘cubic’

雙三次插值'v4'-MATLAB提供的插值方法缺省時,雙線性插值例8在某海域測得一些點(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）×（-50，150）里的哪些地方船要避免進入．4.作出水深小于5的海域范圍,即z=5的等高線...3作海底曲面圖x=xlsread('hd.xls','sheet2','A1:N1');y=xlsread('hd.xls','sheet2','A2:N2');z=xlsread('hd.xls','sheet2','A3:N3');z=-z;xi=75:0.5:200;yi=-50:0.2:150;zi=griddata(x,y,z,xi',yi,'cubic')mesh(xi,yi,zi)figurecontour(xi,yi,zi,[-5-5])gridonholdonplot(x,y,'*')

數(shù)據(jù)擬合

擬合問題引例1溫度t(oC)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求60oC時的電阻R．

設

R=at+ba,b為待定系數(shù)曲線擬合問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個點（xi,yi)i=1,…,n,

尋求一個函數(shù)（曲線）y=f(x),

使f(x)

在某種準則下與所有數(shù)據(jù)點最為接近，即曲線擬合得最好．

+++++++++xyy=f(x)(xi,yi)

i

i為點（xi,yi)與曲線y=f(x)的距離擬合與插值的關系

函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構造一個函數(shù)作為近似，由于近似的要求不同，二者在數(shù)學方法上是完全不同的．實例1：下面數(shù)據(jù)是某次實驗所得，希望得到X和f之間的關系？問題：給定一批數(shù)據(jù)點，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點，而是要求它反映對象整體的變化趨勢，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合．若要求所求曲線（面）通過所給所有數(shù)據(jù)點，就是插值問題；最臨近插值、線性插值、樣條插值與曲線擬合結果：曲線擬合問題最常用的解法——線性最小二乘法的基本思路第一步:先選定一組函數(shù)

r1(x),r2(x),…,rm(x),m<n,

令

f(x)=a1r1(x)+a2r2(x)+…+amrm(x)（1）其中

a1,a2,…,am

為待定系數(shù)．

第二步:確定a1,a2,…,am

的準則（最小二乘準則）：使n個點（xi,yi)與曲線y=f(x)的距離

i

的平方和最小

．記

問題歸結為，求

a1,a2,…,am

使

J(a1,a2,…,am)

最?。€性最小二乘擬合f(x)=a1r1(x)+…+amrm(x)中函數(shù){r1(x),…,rm(x)}的選取

1.通過機理分析建立數(shù)學模型來確定f(x)；

2.將數(shù)據(jù)(xi,yi)i=1,…,n

作圖，通過直觀判斷確定f(x)：++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx用MATLAB解擬合問題1.線性最小二乘擬合2.非線性最小二乘擬合用MATLAB作線性最小二乘擬合1.作多項式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.多項式在x處的值y可用以下命令計算：

y=polyval（a，x）輸出擬合多項式系數(shù)a=[a1,…,am,

am+1](數(shù)組)）輸入同長度的數(shù)組x，y擬合多項式次數(shù)即要求出二次多項式:中的使得:例2對下面一組數(shù)據(jù)作二次多項式擬合2）計算結果：Ａ=-9.810820.1293-0.0317x=0.1:0.1:1.1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];a=polyfit(x,y,2);y1=polyval(a,x);plot(x,y,'o',x,y1)holdona=polyfit(x,y,3);y1=polyval(a,x);plot(x,y1,'r--')1.lsqcurvefit已知數(shù)據(jù)點：

xdata=（xdata1,xdata2,…,xdatan）,ydata=（ydata1,ydata2,…,ydatan）

用MATLAB作非線性最小二乘擬合

MATLAB提供了兩個求非線性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin．兩個命令都要先建立M文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x,xdata1）,…,F（x，xdatan））T中的參變量x(向量),使得

輸入格式為:

(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);

(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,

ydata,…);fun是一個事先建立的定義函數(shù)F(x,xdata)

的M文件,自變量為x和xdata迭代初值已知數(shù)據(jù)點選項見無約束優(yōu)化說明:x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得