05插值和擬合
黃河小浪底調水調沙問題
data3.txt
1800 1900 2100 2200 2300 2400 2500 2600 2650 2700 2720 2650 32 60 75 85 90 98 100 102 108 112 115 116 2600 2500 2300 2200 2000 1850 1820 1800 1750 1500 1000 900 118 120 118 105 80 60 50 30 26 20 8 5
插值
% 問題一 % v:水流量S:含沙量；V:排沙量? % 假設水流量和含沙量都是連續的，某一時刻的排沙量V=v(t)S(t) % 先已知某些時刻的 水流量 和 含沙量 ，給出估計任意時刻的 排沙量 及 總排沙量 % 總排沙量 是對 排沙量 做 積分 % 時間8：00-20：00 t1=28800,t24=1022400=t2 format compact; clc,clear; load data3.txt liu = data3([1,3],:); liu=liu';liu=liu(:); % 提出水流量並按照順序變成列向量 sha = data3([2,4],:); sha=sha';sha=sha(:); % 提出含沙量並按照順序變成列向量 y=sha.*liu;y=y'; % 計算排沙量，變成行向量 i=1:24; t=(12*i-4)*3600; t1=t(1);t2=t(end); % 插值 pp=csape(t,y); % 進行三次樣條插值 xsh=pp.coefs % 求得插值多項式的系數矩陣,每一行是一個區間上多項式的系數 TL = quadl(@(tt)fnval(pp,tt),t1,t2) % 求 總排沙量 的積分運算 % 可視化展示插值多項式和原來的值 t0=t1:0.1:t2; y0=fnval(pp,t0); plot(t,y,'+',t0,y0)
xsh =
1.0e+05 *
-0.0000 -0.0000 0.0000 0.5760
-0.0000 -0.0000 0.0000 1.1400
-0.0000 -0.0000 0.0000 1.5750
0.0000 -0.0000 0.0000 1.8700
-0.0000 0.0000 0.0000 2.0700
0.0000 -0.0000 0.0000 2.3520
0.0000 0.0000 0.0000 2.5000
-0.0000 0.0000 0.0000 2.6520
0.0000 -0.0000 0.0000 2.8620
-0.0000 0.0000 0.0000 3.0240
0.0000 -0.0000 0.0000 3.1280
-0.0000 0.0000 -0.0000 3.0740
-0.0000 -0.0000 0.0000 3.0680
0.0000 -0.0000 -0.0000 3.0000
-0.0000 0.0000 -0.0000 2.7140
0.0000 -0.0000 -0.0000 2.3100
0.0000 0.0000 -0.0000 1.6000
-0.0000 0.0000 -0.0000 1.1100
0.0000 -0.0000 -0.0000 0.9100
-0.0000 0.0000 -0.0000 0.5400
0.0000 -0.0000 -0.0000 0.4550
0.0000 -0.0000 -0.0000 0.3000
0.0000 0.0000 -0.0000 0.0800
TL =
1.8440e+11
擬合
% 問題二：確定排沙量和水流量的關系
format compact;
% 畫出排沙量和水流量的散點圖
clc,clear;
load data3.txt
liu = data3([1,3],:); liu=liu';liu=liu(:); % 提出水流量並按照順序變成列向量
sha = data3([2,4],:); sha=sha';sha=sha(:); % 提出含沙量並按照順序變成列向量
y=sha.*liu; % 計算排沙量，這里是列向量
subplot(1,2,1),plot(liu(1:11),y(1:11),'*')
subplot(1,2,2),plot(liu(12:24),y(12:24),'*')
% 第一階段基本上是線性關系
% 第一階段和第二階段都准備用一次和二次曲線擬合
% 哪個模型的剩余標准差小就選取哪個模型
format long e
% 以下是第一階段的擬合
for j=1:2
nihe1{j}=polyfit(liu(1:11),y(1:11),j); % 擬合多項式,系數排列從高次冪到低次冪
yhat1{j}=polyval(nihe1{j},liu(1:11)); % 求預測值
cha1(j)=sum((y(1:11)-yhat1{j}).^2); % 求誤差平方和
rmse1(j)=sqrt(cha1(j)/(10-j)); % 求剩余標准差
end
celldisp(nihe1) % 顯示細胞數組的所有元素
rmse1
% 以下是第二階段的擬合
for j=1:2
nihe2{j}=polyfit(liu(12:24),y(12:24),(j)); % 使用細胞數組
yhat2{j}=polyval(nihe2{j},liu(12:24)); % 求預測值
cha2(j)=sum((y(12:24)-yhat2{j}).^2); % 求誤差平方和
rmse2(j)=sqrt(cha2(j)/(11-j)); % 求剩余標准差
end
celldisp(nihe2) % 顯示細胞數組的所有元素
rmse2
format % 恢復默認短小數的顯示格式
% 最終結果
% 一: y=250.5655v-373384.4661
% 二: y=0.167v*2-180.4668v+72421.0982
