数学建模(一):插值

目录

题目:

实现代码:

Tips:


题目:

实现代码:

%%
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)的积分

 

本图文内容来源于网友网络收集整理提供,作为学习参考使用,版权属于原作者。
THE END
分享
二维码
< <上一篇
下一篇>>