數學模型-插值與拟合matlab實作

文章目錄

• • • • 插值與拟合模型(一)----水道測量問題
      • • 問題1 插值實驗
        • 問題2 水道測量問題
      • 插值與拟合模型(二)-----水塔流量問題

插值與拟合模型(一)----水道測量問題

Matlab的一維插值函數為interp1()， 調用格式為：

yy=interp1(x,y,xx,方法) 其中x=[x1,x2,…,xn]’, y=[y1,y2,…,yn]’,

兩個向量分别為一組自變量和函數值

xx為待求插值點處橫坐标， yy傳回的對應縱坐标

插值方法可選用方式
’ linear’	線性插值
’ nearest’	最近鄰等值方式
'cubic’	三次Hermite插值
’ spline’	三次樣條插值

問題1 插值實驗

Matlab程式chazhi1.m：

x=0:0.1:1;
y=(x.^2-3*x+7).*exp(-4*x).*sin(2*x); %産生原始資料
subplot(1,2,1);%一行兩個圖，此處為第一個圖
plot(x,y,x,y,'ro') %作圖
xx=0:0.02:1; %待求插值點
yy=interp1(x,y,xx,'spline'); %此處可用nearest,cubic,spline分别試驗
%或yy=spline(x,y,xx);
subplot(1,2,2)%一行兩個圖，此處為第二個圖
plot(x,y,'ro',xx,yy,'b') %作圖
           

問題2 水道測量問題

解答：

所給14個點平面散點圖， 其中有兩點不在區域： (75,200)ⅹ (-50,150)

反距離權重法 (IDW)算法

設有 n n n 個點 ( x i , y i , z i ) \left(x_{i}, y_{i}, z_{i}\right) (xi​,yi​,zi​), 計算平面上任意點 ( x , y ) (x, y) (x,y) 的 z z z 值。

z = ∑ i = 1 n w i ⋅ z i z=\sum_{i=1}^{n} w_{i} \cdot z_{i} z=i=1∑n​wi​⋅zi​

其中權重: w i = 1 / d i p ∑ i = 1 n 1 / d i p , d i = ( x − x i ) 2 + ( y − y i ) 2 w_{i}=\frac{1 / d_{i}^{p}}{\sum_{i=1}^{n} 1 / d_{i}^{p}}, d_{i}=\sqrt{\left(x-x_{i}\right)^{2}+\left(y-y_{i}\right)^{2}} wi​=∑i=1n​1/dip​1/dip​​,di​=(x−xi​)2+(y−yi​)2

​

當 p p p 越大, 則當 ( x i , y i ) \left(x_{i}, y_{i}\right) (xi​,yi​) 距離 ( x , y ) (x, y) (x,y) 越近, 其相對作用越大，越遠相對作用越小。

這裡取 p = 3 p=3 p=3 。

按照 IDW 方法, 對 x x x 方向和 y y y​​ 方向都按照間隔 1 進行剖分。 得到的水底河床圖及等值線圖

Matlab程式如下：

clear; %AMCM86A
data=[129.0,7.5,4;
140.0,141.5,8;
108.5,28.0,6;
88.0,147.0,8;
185.5,22.5,6;
195.0,137.5,8;
105.5,85.5,8;
157.5,-6.5,9;
107.5,-81.0,9;
77.0,3.0,8;
81.0,56.5,8;
162.0,84.0,4;
117.5,-38.5,9;
162.0,-66.5,9];%水道測量資料
plot(data(:,1),data(:,2),'*')
xx=[75,-50;
200,-50;
200,150;
75,150;
75,-50];
hold on;
plot(xx(:,1),xx(:,2),'r')
xlabel('X');
ylabel('Y');
title('水道離散點平面圖');
pause
d=zeros(14,1);
w=zeros(14,1);
p=3; %參數p
[X,Y]=meshgrid(75:1:200,-50:1:150);
%對(X,Y)進行網格剖分
Z=zeros(size(X));
[m,n]=size(X);
tp=0;
for i=1:m
for j=1:n
for k=1:14
%(x,y)格點到各已知各點距離
d(k)=sqrt((X(i,j)-data(k,1))^2+(Y(i,j)-data(k,2))^2);
w(k)=1.0/d(k)^p;%各點權重
end;
s=sum(w);
w=w/s;%權值歸一化
z=sum(data(:,3).*w);
%權重求和
Z(i,j)=-z;
if z<=5
tp=tp+1;
D(tp,1)=X(i,j); D(tp,2)=Y(i,j);
%獲得水深低于5英尺的點
end
end
end
subplot(2,1,1);
surf(X,Y,Z) %作曲面圖
xlabel('X'); ylabel('Y'); zlabel('Z');
title('水道測量水底形狀圖');
subplot(2,1,2);
contour(X,Y,Z); %作等值線圖
hold on
plot(D(:,1),D(:,2),'*')%作出水深低于5英尺的點
xlabel('X'); ylabel('Y');
title('水道測量等值線圖');
grid on
hold off
           

插值與拟合模型(二)-----水塔流量問題

（ MCM91A） 美國某洲的各用水管理機構要求各社群提供以每小時多少加侖計的用水率，以及每天總的用水量，但許多社群并沒有測量水流入或流出水塔水量的裝置，他們隻能每小時測量水塔中的水位，精度在0.5%以内，更為重要的是，無論什麼時候，隻要水塔中的水位下降到某一最低水位L時，水泵就啟動向水塔重新充水至某一最高水位H，但也無法得到水泵的供水量的測量資料。 水泵每天向水塔充水一次或兩次， 每次約兩小時。水塔是一個垂直圓形柱體， 高為40英尺， 直徑57英尺。當水塔的水位降至27.00英尺時開始向水塔充水，當水位升至35.50英尺時停止充水。

試估計在任何時刻, 甚至包括水泵正在工作期間内, 水從水塔流出的流量 f ( t ) f(t) f(t), 并估計一天的總用水量, 表 1 中給出了某個真實小鎮某一天的真實資料。

解答:

1. 水塔充水時間的确定

   (1) 第一次充水時間： t=32284到t=39435秒， 間隔dt=1.9864小時

   (2) 第二次充水時間： t=75021到t=82649秒， 間隔dt=2.1189小時

2. 計算各時刻塔内水的體積

   機關轉換為1英尺=0.3048米, 1升=1/3.785411加侖

   體積計算公式為 v = π ⋅ d 2 h / 4 v=\pi \cdot d^{2} h / 4 v=π⋅d2h/4​

3.計算各時刻點的水流量(加侖/小時)

水流量公式為: f ( t ) = ∣ d v ( t ) d t ∣ \quad f(t)=\left|\frac{d v(t)}{d t}\right| f(t)=∣∣∣​dtdv(t)​∣∣∣​​

25 個時刻處的水流量采用差分的方法得到，共分三段分别處理

差分公式為：

(1) 對每段前兩點采用向前三點差分公式

f ( t i ) = ∣ − 3 V i + 4 V i + 1 − V i + 2 2 ( t i + 1 − t i ) ∣ f\left(t_{i}\right)=\left|\frac{-3 V_{i}+4 V_{i+1}-V_{i+2}}{2\left(t_{i+1}-t_{i}\right)}\right| f(ti​)=∣∣∣∣​2(ti+1​−ti​)−3Vi​+4Vi+1​−Vi+2​​∣∣∣∣​

(2) 對每段最後兩點采用向後三點差分公式

f ( t i ) = ∣ 3 V i − 4 V i − 1 + V i − 2 2 ( t i − t i − 1 ) ∣ f\left(t_{i}\right)=\left|\frac{3 V_{i}-4 V_{i-1}+V_{i-2}}{2\left(t_{i}-t_{i-1}\right)}\right| f(ti​)=∣∣∣∣​2(ti​−ti−1​)3Vi​−4Vi−1​+Vi−2​​∣∣∣∣​

(3) 對每段中間點采用中心差分公式

f ( t i ) = ∣ − V i + 2 + 8 V i + 1 − 8 V i − 1 + V i − 2 12 ( t i + 1 − t i ) ∣ \quad f\left(t_{i}\right)=\left|\frac{-V_{i+2}+8 V_{i+1}-8 V_{i-1}+V_{i-2}}{12\left(t_{i+1}-t_{i}\right)}\right| f(ti​)=∣∣∣∣​12(ti+1​−ti​)−Vi+2​+8Vi+1​−8Vi−1​+Vi−2​​∣∣∣∣​

計算各點水流量見表3.

4.用三次樣條拟合流量資料

對表3中25個時刻點流量數

據采用三次樣條插值見圖6

5.**一天總用水量計算 **

方法 1:直接積分法:

S 1 = ∫ 0 24 f ( t ) d t = 332986 ( S_{1}=\int_{0}^{24} f(t) d t=332986\left(\right. S1​=∫024​f(t)dt=332986( 加侖） \quad

\quad ​​ 設樣條揷值得 n \mathrm{n} n​​ 個點 y 1 , y 2 , … , y n y_{1}, y_{2}, \ldots, y_{n} y1​,y2​,…,yn​​​, 間隔為 h h h​​ 數值積分公式: S = [ y 1 + y n + 2 ( y 2 + … + y n − 1 ) ] ⋅ h 2 S=\left[y_{1}+y_{n}+2\left(y_{2}+\ldots+y_{n-1}\right)\right] \cdot \frac{h}{2} S=[y1​+yn​+2(y2​+…+yn−1​)]⋅2h​​

方法2：分段計算

第一次充水前用水 V 1 = 606125 − 514872 = 91253 V_{1}=606125-514872= 91253 V1​=606125−514872=91253​​​​ (加侖)

第二次充水前用水 V 2 = 677715 − 514872 = 162843 V_{2}=677715-514872=162843 V2​=677715−514872=162843​​( 加侖 ) ) )​

[ 22.9581 , 23.88 ] [22.9581,23.88] [22.9581,23.88]​​ 期間用水 V 3 = 677715 − 663397 = 14318 V_{3}=677715-663397=14318 V3​=677715−663397=14318​​ (加侖 )

第一次充水期間用水: V 4 = ∫ 8.9678 10.9542 f ( t ) d t = 30326 V_{4}=\int_{8.9678}^{10.9542} f(t) d t=30326 V4​=∫8.967810.9542​f(t)dt=30326​ (加侖)

第二次充水期間用水: V 5 = ∫ 20.8392 22.9581 f ( t ) d t = 31605 V_{5}=\int_{20.8392}^{22.9581} f(t) d t=31605 V5​=∫20.839222.9581​f(t)dt=31605​​( 加侖 )

[ 23.88 , 24 ] [23.88, 24] [23.88,24]期間用水: V 6 = ∫ 23.88 24 f ( t ) d t = 1524 V_{6}=\int_{23.88}^{24} f(t) d t=1524 V6​=∫23.8824​f(t)dt=1524​​​​​​​​ (加侖)

總共用水 S 2 = ∑ i = 1 6 V i = 33186 S_{2}=\sum_{i=1}^{6} V_{i}=33186 S2​=∑i=16​Vi​=33186 (加侖)

兩種方法誤差 err ⁡ = ∣ S 1 − S 2 S 1 ∣ = 0.34 % \operatorname{err}=\left|\frac{S_{1}-S_{2}}{S_{1}}\right|=0.34 \% err=∣∣∣​S1​S1​−S2​​∣∣∣​=0.34%

6.水泵流量計算

第一次充水期間水塔體積增加	充水時間	第一次充水期間水泵平均流量
Δ V 1 = 677715 − 514872 = 162843 \Delta V_{1}=677715-514872=162843 ΔV1​=677715−514872=162843(加侖)	Δ t 1 = 10.9542 − 8.9678 = 1.9864 \Delta t_{1}=10.9542-8.9678=1.9864 Δt1​=10.9542−8.9678=1.9864( 小時 )	$ p_{1}=\frac{\Delta V_{1}+\int_{8.9678}^{10.9542} f(t) d t}{\Delta t_{1}}=97246$(加侖/小時)

第二次充水期間水塔體積增加	充水時間	第二次充水期間水泵平均流量
Δ V 2 = 677715 − 514872 = 162843 \Delta V_{2}=677715-514872=162843 ΔV2​=677715−514872=162843(加侖)	\quad Δ t 2 = 22.9581 − 20.8392 = 2.1189 \Delta t_{2}=22.9581-20.8392=2.1189 Δt2​=22.9581−20.8392=2.1189(小時 )	p 2 = Δ V 2 + ∫ 20.8392 22.9581 f ( t ) d t Δ t 2 = 91769 p_{2}=\frac{\Delta V_{2}+\int_{20.8392}^{22.9581} f(t) d t}{\Delta t_{2}}=91769 p2​=Δt2​ΔV2​+∫20.839222.9581​f(t)dt​=91769( 加侖 小時 )

則整個充水期間水泵平均流量

p = ( p 1 + p 2 ) / 2 = 94507 p=\left(p_{1}+p_{2}\right) / 2=94507 p=(p1​+p2​)/2=94507​( 加侖 /小時 )

Matlab程式

%AMCM91A
c=0.3048; %1英尺等于0.3048米
p=1.0/3.785; %1升=1/3.785411加侖
d=57*c; %直徑
h=31.75*c;
v=pi*d*d*h/4*1000*p;
data=[0,3175;
3316,3110;
6635,3054;
10619,2994;
13937,2947;
17921,2892;
21240,2850;
25223,2797;
28543,2752;
32284,2697;
39435,3550;
43318,3445;
46636,3350;
49953,3260;
53936,3167;
57254,3087;
60574,3012;
64554,2927;
68535,2842;
71854,2767;
75021,2697;
82649,3550;
85968,3475;
89953,3397;
93270,3340];%原始資料
t=data(:,1)/3600;%計算時間(小時為機關)
v=pi*d*d*data(:,2)/100*c/4*1000*p;%計算體積
%計算差分
n=length(v);
f=zeros(n,1); %存儲差分值
%計算第一段
n1=10;
for i=1:n1
if i<=2 %前兩點采用向前三點差分
f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));
elseif i<=n1-2
%采用中心差分公式
f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));
elseif i>=n1-1 %後兩點采用向後三點差分
f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));
end
end
n2=21; %計算第二段
for i=n1+1:n2
if i<=n1+2 %前兩點采用向前差分
f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));
elseif i<=n2-2
f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));
elseif i>=n2-1
f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));
end
end
n3=25; %計算第三段
for i=n2+1:n3
if i<=n2+2 %前兩點采用向前差分
f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));
elseif i<=n3-2
f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));
elseif i>=n3-1
f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));
end
end
plot(t,f,'r*'); %畫原始點圖
tmin=min(t); tmax=max(t);
tt=tmin:0.1:tmax; %獲得離散的時間點,用于作樣條曲線
ff=spline(t,f,tt); %計算三次樣條插值
hold on
plot(tt,ff,'b'); %畫樣條曲線
xlabel('時間(小時)');
ylabel('流量(加侖/小時)');
title('水塔流量圖');
hold off
dt=0.05;
t2=0.5:dt:24.5; %獲得離散的時間點,用于積分
nn=length(t2);
f2=spline(t,f,t2);
%計算24小時用水量,采用複化梯形公式
s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)))*dt/2 ;
fprintf('(全部積分法)1天總水流量s= %8.2f\n',s);
%第一次水塔增加的水
v10=v(11)-v(10);
dt1=t(11)-t(10);
%第一次充水其間流出的水
tp=t(10):dt:t(11);
nn=length(tp);
yp=spline(t,f,tp); %計算三次樣條插值
v11=(yp(1)+yp(nn)+2*sum(yp(2:nn-1)))*dt/2;
v1=v10+v11;%第一次充水總量
p1=v1/dt1; %第一次充水的平均水流量
%第二次水塔增加的水
v20=v(22)-v(21);
dt2=t(22)-t(21);
%第二次充水其間流出的水
tp1=t(21):dt:t(22);
nn=length(tp1);
yp1=spline(t,f,tp1); %計算三次樣條插值
v21=(yp1(1)+yp1(nn)+2*sum(yp1(2:nn-1)))*d
v2=v20+v21;%第二次充水總量
p2=v2/dt2; %第二次充水的平均水流量
p=(p1+p2)/2; %兩次充水準均水流量
fprintf('兩次充水準均水流量p= %8.2f\n',p);
%第一次充水前總流量
vv1=v(1)-v(10);
%兩次充水間總流量
vv2=v(11)-v(21);
% t=83649到85968其間流量
vv3=v(22)-v(23);
%第一次充水期間流量
ta=t(10):dt:t(11);
%獲得離散的時間點,用于積分
nn=length(ta);
fa=spline(t,f,ta);
s1=(fa(1)+fa(nn)+2*sum(fa(2:nn-1)))*dt/2;
%第二次充水期間流量
tb=t(21):dt:t(22); %獲得離散的時間點,用于積分
nn=length(tb);
fb=spline(t,f,tb);
s2=(fb(1)+fb(nn)+2*sum(fb(2:nn-1)))*dt/2;
%t=85968到86400流量
tc=t(23):dt:24; %獲得離散的時間點,用于積分
nn=length(tc);
fc=spline(t,f,tc);
s3=(fc(1)+fc(nn)+2*sum(fc(2:nn-1)))*dt/2 ;
ss=vv1+vv2+vv3+s1+s2+s3;
fprintf('(部分積分法)1天總水流量ss= %8.2f\n',ss);
err=abs((s-ss)/s);
fprintf('兩種計算法總量相對誤差%6.2f%%\n',err*100);
           

程式計算結果為：

(全部積分法)1天總水流量s= 332986.22

兩次充水準均水流量 p= 94507.48

(部分積分法)1天總水流量 ss= 331869.29

兩種計算法總量相對誤差 0.34%