定義
插值和擬合：
曲線擬合是指您擁有散點數據集並找到最適合數據一般形狀的線（或曲線）。
插值是指您有兩個數據點並想知道兩者之間的值是什么。中間的一半是他們的平均值，但如果你只想知道兩者之間的四分之一，你必須插值。
擬合
我們着手寫一個線性方程圖的擬合：
y=3x^3+2x^2+x+2
首先我們生成一組數據來分析：
x=-5:0.5:5; e=50*rand(1,length(x))-25;%制造[-25,25]的隨機數作為誤差 y=3*x.^3+2*x.^2+x+2+e;%得到y值 plot(x,y,'.')
x =
Columns 1 through 6
-5.0000 -4.5000 -4.0000 -3.5000 -3.0000 -2.5000
Columns 7 through 12
-2.0000 -1.5000 -1.0000 -0.5000 0 0.5000
Columns 13 through 18
1.0000 1.5000 2.0000 2.5000 3.0000 3.5000
Columns 19 through 21
4.0000 4.5000 5.0000
y =
Columns 1 through 6
-350.0110 -248.6360 -169.3421 -89.5653 -88.2298 -57.7238
Columns 7 through 12
-32.5505 2.3308 11.5861 9.0123 -0.4538 5.7254
Columns 13 through 18
-2.1840 30.3596 20.4478 73.2138 88.1756 152.0492
Columns 19 through 21
236.2809 334.3864 411.0563
這時候x,y作為已知數據存在，要求我們擬合x,y的散點圖，這時候會用到polyfit這個函數。
語法
p = polyfit(x,y,n)
[p,S] = polyfit(x,y,n)
[p,S,mu] = polyfit(x,y,n)
說明
p = polyfit(x,y,n) 返回階數為 n 的多項式 p(x) 的系數，該階數是 y 中數據的最佳擬合（在最小二乘方式中）。p 中的系數按降冪排列，p 的長度為 n+1
[p,S] = polyfit(x,y,n) 還返回一個結構體 S，后者可用作 polyval 的輸入來獲取誤差估計值。
[p,S,mu] = polyfit(x,y,n) 還返回 mu，后者是一個二元素向量，包含中心化值和縮放值。mu(1) 是 mean(x)，mu(2) 是 std(x)。使用這些值時，polyfit 將 x 的中心置於零值處並縮放為具有單位標准差
plot(x,y,'.')%已知數據的散點圖 hold on d=polyfit(x,y,2);%進行線性擬合，假設函數為二次函數 y1=d(1)*x.^2+d(2)*x+d(3); plot(x,y1)//繪制猜測的函數
明顯不符合，我們繼續猜測
plot(x,y,'.') hold on d=polyfit(x,y,3); y1=d(1)*x.^3+d(2)*x.^2+d(3)*x+d(4); plot(x,y1)
擬合成功。
插值
一維插值
以上面的數據為例子
x=-5:0.5:5; e=50*rand(1,length(x))-25;%制造[-25,25]的隨機數作為誤差 y=3*x.^3+2*x.^2+x+2+e;%得到y值 %插值 x2=-5:0.1:5; y2=interp1(x,y,'x2,'spline'); plot(x2,y2,'r+')
二維插值
語法
[X,Y] = meshgrid(x,y)
[X,Y] = meshgrid(x)
[X,Y,Z] = meshgrid(x,y,z)
[X,Y,Z] = meshgrid(x)
說明
[X,Y] = meshgrid(x,y) 基於向量 x 和 y 中包含的坐標返回二維網格坐標。X 是一個矩陣，每一行是 x 的一個副本；Y 也是一個矩陣，每一列是 y 的一個副本。坐標 X 和 Y 表示的網格有 length(y) 個行和 length(x) 個列。
[X,Y] = meshgrid(x) 與 [X,Y] = meshgrid(x,x) 相同，並返回網格大小為 length(x)×length(x) 的方形網格坐標。
[X,Y,Z] = meshgrid(x,y,z) 返回由向量 x、y 和 z 定義的三維網格坐標。X、Y 和 Z 表示的網格的大小為 length(y)×length(x)×length(z)。
[X,Y,Z] = meshgrid(x) 與 [X,Y,Z] = meshgrid(x,x,x) 相同，並返回網格大小為 length(x)×length(x)×length(x) 的三維網格坐標。
我們先繪制一個三維的圖像
[x,y]=meshgrid(-5:1:5,-10:2:10); z=12-x.^3-y.^3; surf(x,y,z)
給人的感覺是很“粗糙”的，這時候我們使用插值，增加數據。
[x0,y0]=meshgrid(-5:0.2:5,-10:0.4:10); z0=interp2(x,y,z,x0,y0,'spline'); surf(x0,y0,z0)
例題
例1
1、 已知以下實驗測量數據，
|
|
1 |
1.5 |
2.5 |
3 |
4.5 |
|
|
-0.86 |
13.66 |
90.25 |
152.57 |
473.70 |
（1） 請畫出以上數據的散點圖；
（2） 請用合適的多項式擬合上述數據；
（3） 請用三次樣條一次插值以上數據.
2、 已知網格化數據如下：
|
|
-1 |
-0.5 |
0 |
0.5 |
1 |
| -2 |
1.77 |
1.12 |
0.32 |
0.44 |
-2.35 |
| -1 |
1.62 |
1.29 |
0.94 |
-1.49 |
-0.78 |
| 0 |
0.23 |
-1.01 |
-0.74 |
1.08 |
-0.13 |
| 1 |
-0.61 |
-0.41 |
-0.63 |
-0.05 |
1.44 |
| 2 |
-2.94 |
-0.22 |
0.56 |
0.17 |
1.73 |
請用三次樣條二次插值法插出更光滑精細的曲面.
答案
%1-(1) clc,clear x=[1 1.5 2.5 3 4.5]; y=[-0.86 13.66 90.25 152.57 473.70]; plot(x,y,'.'); %1-(2) d=polyfit(x,y,1); y1=d(1)*x+d(2); plot(x,y,'.'); plot(x,y1);%對比線性圖與散點圖，發現不匹配 d=polyfit(x,y,2); y1=d(2)*x.^2+d(2)*x+d(3); plot(x,y,'.'); plot(x,y1);%對比線性圖與散點圖，發現匹配成功，擬合完成 %1-(3) x=[1 1.5 2.5 3 4.5]; y=[-0.86 13.66 90.25 152.57 473.70]; x1=[1:0.1:4.5]; y1=interp1(x,y,x1,'spline'); plot(x,y,x1,y1); %2 [x1,y1]=meshgrid(-2:1:2,-1:0.5:1); z1=[1.77 1.21 0.32 0.44 -2.35;1.62 1.29 0.94 -1.49 -0.78;0.23 -1.01 -0.74 1.08 -0.13;... -0.61 -0.41 -0.63 -0.05 1.44;-2.94 -0.22 0.56 0.17 1.73]; surf(x1,y1,z1) [x,y]=meshgrid(-1:0.1:1, -2:0.2:2); z=interp2(x1,y1,z1,x,y,'spline'); surf(x,y,z)
例2
數據：點擊下載
實際任務就是sheet2中中的元素，只給了單周的元素數據，沒有給偶數周，這就需要我們插值出來
答案
clc,clear
%% 讀取數據
ydatas = xlsread('data1.xls',2,'C5:J17');
% 清除無效數據
ydatas(9:10,:)=[];
%% 繪圖准備
% y坐標軸名稱
yname = ["周數","輪蟲(10^6/L)","溶氧(mg/l)","COD(mg/l)","水溫(℃)","PH值","鹽度","透明度(cm)","總鹼度","氯離子","透明度","生物量"];
% yname = {'輪蟲(10^6/L)','溶氧(mg/l)','COD(mg/l)','水溫(℃)','PH值','鹽度','透明度(cm)','總鹼度','氯離子','透明度','生物量'};
% x坐標
xdatas = 1:2:15;
xval = 1:15;
%% 線性插值
figure(1)
yval1 = [];
for index=1:11
tmp = interp1(xdatas,ydatas(index,:),xval,'linear');
yval1 = [yval1;tmp];
subplot(3,4,index)
plot(xdatas,ydatas(index,:),'*',xval,tmp,'--r')
xlabel('時間/周');
ylabel(yname(index+1))
legend('原始數據','線性插值數據')
end
%% 三次樣條插值
figure(2)
yval2 = [];
for index=1:11
tmp = interp1(xdatas,ydatas(index,:),xval,'spline');
yval2 = [yval2;tmp];
subplot(3,4,index)
plot(xdatas,ydatas(index,:),'*',xval,tmp,'--r')
xlabel('時間/周');
ylabel(yname(index+1))
legend('原始數據','三次樣條插值數據')
end
%% 三次多項式插值
figure(3)
yval3 = [];
for index=1:11
tmp = interp1(xdatas,ydatas(index,:),xval,'pchip');
yval3 = [yval3;tmp];
subplot(3,4,index)
plot(xdatas,ydatas(index,:),'*',xval,tmp,'--r')
xlabel('時間/周');
ylabel(yname(index+1))
legend('原始數據','三次多項式插值數據')
end
%% 數據存儲
new_data1 = [xval;yval1];
new_data1 = [yname',new_data1];
xlswrite('output.xls',new_data1,'Sheet1');
new_data2 = [xval;yval2];
new_data2 = [yname',new_data2];
xlswrite('output.xls',new_data2,'Sheet2');
new_data3 = [xval;yval3];
new_data3 = [yname',new_data3];
xlswrite('output.xls',new_data3,'Sheet3');
輸出結果數據：https://www.lanzous.com/iau54la
例3
數據下載：點擊下載
答案
clc,clear
%%
xdatas = xlsread('data2.xlsx','A2:A11');
ydatas = xlsread('data2.xlsx','B2:B11');
a=polyfit(xdatas,ydatas,1);
y1=a(1)*xdatas+a(2);
subplot(3,2,1)
plot(xdatas,y1,'r-',xdatas,ydatas,'*')
title('一次擬合');
xlabel('年份');
ylabel('人口/萬');
y1_datas = a(1)*xdatas+a(2);
ssr1 = sum((y1_datas-mean(ydatas)).^2);
sst1 = sum((ydatas-mean(ydatas)).^2);
r1_sq = ssr1/sst1;
txt = ['R^2=' num2str(r1_sq)];
text(2010,138000,txt)
%%
b=polyfit(xdatas,ydatas,2);
y2=b(1)*xdatas.^2+b(2)*xdatas+b(3);
subplot(3,2,2)
plot(xdatas,y2,'r-',xdatas,ydatas,'*')
title('二次擬合');
xlabel('年份');
ylabel('人口/萬');
y2_datas = b(1)*xdatas.^2+b(2)*xdatas+b(3);
ssr2 = sum((y2_datas-mean(ydatas)).^2);
sst2 = sum((ydatas-mean(ydatas)).^2);
r2_sq = ssr2/sst2;
txt = ['R^2=' num2str(r2_sq)];
text(2010,138000,txt)
%%
c=polyfit(xdatas,ydatas,3);
y3=c(1)*xdatas.^3+c(2)*xdatas.^2+c(3)*xdatas+c(4);
subplot(3,2,3)
plot(xdatas,y3,'r-',xdatas,ydatas,'*')
title('三次擬合');
xlabel('年份');
ylabel('人口/萬');
y3_datas = c(1)*xdatas.^3+c(2)*xdatas.^2+c(3)*xdatas+c(4);
ssr3 = sum((y3_datas-mean(ydatas)).^2);
sst3 = sum((ydatas-mean(ydatas)).^2);
r3_sq = ssr3/sst3;
txt = ['R^2=' num2str(r3_sq)];
text(2010,138000,txt)
%%
d=polyfit(xdatas,ydatas,4);
y4=d(1)*xdatas.^4+d(2)*xdatas.^3+d(3)*xdatas.^2+d(4)*xdatas+d(5);
subplot(3,2,4)
plot(xdatas,y4,'r-',xdatas,ydatas,'*')
title('四次擬合');
xlabel('年份');
ylabel('人口/萬');
y4_datas = d(1)*xdatas.^4+d(2)*xdatas.^3+d(3)*xdatas.^2+d(4)*xdatas+d(5);
ssr4 = sum((y4_datas-mean(ydatas)).^2);
sst4 = sum((ydatas-mean(ydatas)).^2);
r4_sq = ssr4/sst4;
txt = ['R^2=' num2str(r4_sq)];
text(2010,138000,txt)
%%
e=polyfit(xdatas,ydatas,5);
y5=e(1)*xdatas.^5+e(2)*xdatas.^4+e(3)*xdatas.^3+e(4)*xdatas.^2+e(5)*xdatas+e(6);
subplot(3,2,5)
plot(xdatas,y5,'r-',xdatas,ydatas,'*')
title('五次擬合');
xlabel('年份');
ylabel('人口/萬');
y5_datas = e(1)*xdatas.^5+e(2)*xdatas.^4+e(3)*xdatas.^3+e(4)*xdatas.^2+e(5)*xdatas+e(6);
ssr5 = sum((y5_datas-mean(ydatas)).^2);
sst5 = sum((ydatas-mean(ydatas)).^2);
r5_sq = ssr5/sst5;
txt = ['R^2：' num2str(r5_sq)];
text(2010,138000,txt)
通過對數據擬合，對比確定系數，反應擬合效果，確定出最佳擬合函數為四次擬合函數：y = 0.26*x^4-2.06*x^3+6.21*x^2-8.32*x+4.19