Matlab如何求离散点的导数
1.通过差分估算
已知同维度的x和y序列，则可使⽤diff(y)./diff(x)来估算。设x为n维向量，Dx=diff(x) 计算向量x的向前差分，DX(i)=X(i+1)-
X(i)，0<i<n。
例⼀
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...
5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...
0.68 0.30];
x=0:0.04:1.48;
plot(x,y)
dy=diff(y)./diff(x);
dx=0.04:0.04:1.48;
figure
plot(dx,dy)
2.通过梯度估算
已知同维度的x和y序列，则可使⽤gradient(y)./gradient(x)来估算 。梯度⽤的是中⼼点差分，diff()⽤的前后两点差分；所以从区间上看梯度⽤的范围⽐导数⼤⼀倍！所以梯度⽅式精度会更⾼⼀些！但是梯度法的边界可能会出现误差问题。
2.1梯度的含义与使⽤
[Fx,Fy]=gradient(F)，其中Fx为其⽔平⽅向上的梯度，Fy为其垂直⽅向上的梯度，Fx的第⼀列元素为
原矩阵第⼆列与第⼀列元素之差，Fx 的第⼆列元素为原矩阵第三列与第⼀列元素之差除以2，以此类推：Fx(i,j)=(F(i,j+1)-F(i,j-1))/2。最后⼀列则为最后两列之差。同理，可以得到Fy。
例⼆
x=[6,9,3,4,0;5,4,1,2,5;6,7,7,8,0;7,8,9,10,0];
[Fx,Fy]=gradient(x)
运⾏结果：
Fx =
3.0000 -1.5000 -2.5000 -1.5000 -
4.0000
-1.0000 -2.0000 -1.0000 2.0000 3.0000
1.0000 0.5000 0.5000 -3.5000 -8.0000
1.0000 1.0000 1.0000 -4.5000 -10.0000
Fy =
-1.0000 -5.0000 -2.0000 -2.0000 5.0000
0 -1.0000 2.0000 2.0000 0
1.0000
2.0000 4.0000 4.0000 -2.5000
1.0000 1.0000
2.0000 2.0000 0
2.2通过梯度估算
例三：差分与梯度对⽐
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...
5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...
0.68 0.30];
x=0:0.04:1.48;
dy=gradient(y)./gradient(x);
dx=0:0.04:1.48;
dx1=0.04:0.04:1.48;
figure
plot(dx,dy,'r')
dy1=diff(y)./diff(x);
hold on
plot(dx1,dy1,'b')
legend('梯度','差分')
运⾏结果
例四
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...
5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...
0.68 0.30];
x=0:0.04:1.48;
p=polyfit(x,y,10); %使⽤10次多项式拟合
y3=polyval(p,x); %求出预测值
plot(x,y,'b',x,y3,'r')
legend('原始数据','拟合函数')
dy=diff(y)./diff(x);
dx=0.04:0.04:1.48;
% dy3=diff(y3)./diff(x);
p1=polyder(p);
dy3=polyval(p1,dx);
figure
plot(dx,dy,'b',dx,dy3,'r')
legend('差分求导','拟合多项式求导')
运⾏结果
4.3点、4点、5点求导法
例五：使⽤5点求导法[5]
y=[7.86 7.84 7.82 7.77 7.72 7.68 7.61 7.51 7.42 7.33 7.21 7.07 6.94 6.79 6.64 6.48 6.29 6.11 ...
5.92 5.72 5.50 5.27 5.03 4.78 4.53 4.25 3.98 3.69 3.40 3.10 2.78 2.43 2.09 1.77 1.42 1.09 ...
0.68 0.30]';
x=(0:0.04:1.48)';
y1 = fivePoint1Order(y,0.04);
figure
plot(x(2:end,1),y1,'r')
function [value] = fivePoint1Order(Sig,h)
% 使⽤五点法求⼀阶导数
matlab求导% Sig被求导的列向量
% h步长，单位为秒
lengSig = size(Sig,1);
value = zeros(lengSig-1,1);
% 边缘点使⽤差分，参考matlab梯度函数的做法
value(1,1) = (Sig(2,1) - Sig(1,1))/h;
for i = 1:lengSig-4
fivePoints = Sig(i:i+4,1);
f_2 = fivePoints(1,1);
f_1 = fivePoints(2,1);
f1 = fivePoints(4,1);
f2 = fivePoints(5,1);
value(i+1,1) = (f_2 - 8*f_1 + 8*f1 - f2)/(12*h);
end
value(i+2,1) = (Sig(i+3,1) - Sig(i+2,1))/h;
value(i+3,1) = (Sig(i+4,1) - Sig(i+3,1))/h;% 边缘点
end
运⾏结果
参考⽂献
[5]RAFATI FARD M, SHARIAT MOHAYMANY A, SHAHRI M. A new methodology for vehicle trajectory reconstruction based on wavelet analysis [J]. Transportation Research Part C: Emerging T
echnologies, 2017, 74(150-67).
发布评论