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
