matlab学习——05插值和拟合(黄河小浪底调水调沙问题)

原作者: [db:作者] 来自: [db:来源] 收藏邀请

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