目录
数学建模专栏:👉🏻数学建模从0到1👈🏻
题目:
编辑
实现代码:
%% t0=data([1:2:end],:);%从第一行开始,每隔两行取一行数据(得到的是时间数据) t0=t0';%对行向量求转置 t0=t0(:);%全部转成列向量 h0=data([2:2:end],:);%同理,取高度数据 h0=h0'; h0=h0(:); H=17.4; V=pi/4*H*H*h0;%计算各时刻的体积 dv = -1*gradient(V,t0);%对各个时刻进行求导,因为斜率是负的,所以这个地方的速度是负数,需要再乘一个-1转正数 n0 = find(h0 == -1);%找出空格数据 n1=[n0(1)-1:n0(2)+1,n0(3)-1:n0(4)+1]; t=t0; t(n1)=[]; dv1=dv; dv1(n1)=[]; %% %这个地方取四个点行不通,因为边界点不连续,然后不连续的情况下求导误差巨大 n2=[n0(1),n0(2),n0(3),n0(4)]; t3=t0; t3(n2)=[]; dv2=dv; dv2(n2)=[]; %% plot(t,dv1,'b*',t3,dv2,'ro'); plot(t,dv1,'b*') plot(t3,dv2,'ro') %% pp=csape(t,dv1);%进行差值 tt=0:0.1:t(end);%给出差值点 fdv=ppval(pp,tt); plot(tt,fdv,'r*'); I=trapz(tt(1:241),fdv(1:241))%计算24h内总流量的数值积分
Tips:
简单讲一下思路:
(1)首先把题目中的数据以矩阵的形式输入到matlab中,分别将时间和时刻这两组数据存进列向量中,对应下面这段代码:【当然,也可以直接自己手动把两组数据输入到列向量中】
t0=data([1:2:end],:);%从第一行开始,每隔两行取一行数据(得到的是时间数据) t0=t0';%对行向量求转置 t0=t0(:);%全部转成列向量 h0=data([2:2:end],:);%同理,取高度数据 h0=h0'; h0=h0(:);
(2)然后,我们要对各个时刻的点进行求导,求导得到的是速度,然后要删掉四个空点(也就是题目中没有给数值的点),然后还要删掉四个分段的点,因为求导的前提条件就是连续,但是对于分段函数的边界点它不是连续的,题中根据速度散点图可知,可以把图分为三段,即有4个边界点,所以一共去除8个不可用点,代码如下所示:
dv = -1*gradient(V,t0);%对各个时刻进行求导,因为斜率是负的,所以这个地方的速度是负数,需要再乘一个-1转正数 n0 = find(h0 == -1);%找出空格数据 n1=[n0(1)-1:n0(2)+1,n0(3)-1:n0(4)+1];%删除8个不可用点 t=t0; t(n1)=[]; dv1=dv; dv1(n1)=[];
(3)最后对数据进行插值,利用csape函数得到三次样条插值的系数,它的结果是存在一个结构体里的,利用ppval函数利用caspe得到的系数计算出插值,并画出图像代码如下所示:
pp=csape(t,dv1);%进行差值 tt=0:0.1:t(end);%给出差值点 fdv=ppval(pp,tt); plot(tt,fdv,'r*');
图:
编辑
(4)利用matlab中的trapz函数对离散的点进行求积分,代码如下所示:
I=trapz(tt(1:241),fdv(1:241))%计算24h内总流量的数值积分
MATLAB函数trapz(x, y, n), 其中y是x的积分, 使用梯形法则逼近函数y = f(x)的积分