对大量实验测得的数据点用Smoothing spline曲线拟合,然后怎样得出拟合后的曲线的导数曲线?
1条回答 默认 最新
- 技术专家团-Joel 2021-11-16 22:21关注
一个简单的办法是直接用三点差分求导,我给你一个函数
function [dydx,x,y] = threePiontDerivative(x,y) x = x(:); y = y(:); n = length(x); m = length(y); if(n~=m) error('x和y长度不一致') end if(n<3) error('x和y长度不够,长度至少要有3') end dx = diff(x); h = mean(dx); if(max(abs(h-dx))>1e-12) warning('x不是等间距的') if(~all(dx>0)) [x,idx] = sort(x); y = y(idx); end [x,ia] = unique(x); y = y(ia); dydx = nonUniformPoint(x,y); else dydx = zeros(size(y)); dydx(2:end-1) = (y(3:end)-y(1:end-2))/(2*h);%首端点 dydx(1) = (-3*y(1)+4*y(2)-y(3))/(2*h);%中间点 dydx(n) = (y(n-2)-4*y(n-1)+3*y(n))/(2*h);%末端点 end end function dy = nonUniformPoint(x,y) dy = zeros(size(y)); n = length(y); dx = diff(x); h1 = dx(1:end-1); h2 = dx(2:end); h11 = h1.*h1; h22 = h2.*h2; h12 = h1.*h2; a = -h2./(h11+h12); b = (h2-h1)./h12; c = h1./(h22+h12); dy(2:n-1) = a.*y(1:n-2)+b.*y(2:n-1)+c.*y(3:n); dy(1) = -(2*h1(1)+h2(1))/(h11(1)+h12(1))*y(1) +... (h1(1)+h2(1))/h12(1)*y(2) +... -h1(1)/(h22(1)+h12(1))*y(3); dy(n) = (h2(end)/(h11(end)+h12(end)))*y(n-2) +... -(h1(end)+h2(end))/h12(end)*y(n-1) +... (h1(end)+2*h2(end))/(h22(end)+h12(end))*y(n); end
然后调用:
x = 0:pi/3:2*pi; y = sin(x); xp = 0:pi/10:2*pi; yp = spline(x,y,xp); figure(1);clf; plot(x,y,'ro',xp,yp,'b-'); dy = threePiontDerivative(xp,yp); hold on plot(xp,dy,'m--') legend('原数据','Spline后','Spline求导')
效果:
当然我还想到了一种更简单的办法!!!, 一并实现了
x = 0:pi/3:2*pi; y = sin(x); xp = 0:pi/10:2*pi; pp = spline(x,y); yp = ppval(pp, xp); % 分段得到的样条插值值 p = pp.coefs; % 获取分段系数 q = zeros(size(p,1),size(p,2)-1); for i = 1:size(p,1) q(i,:) = polyder(p(i,:)); %直接对分段多项式系数求导 end pp.coefs = q; pp.order = 3; % 把pp更改为求导后的系数和阶数 dy = ppval(pp, xp); % 直接用ppval获取样条插值后的导数 figure(1);clf; % 画图 plot(x,y,'ro',xp,yp,'b-');hold on; plot(xp,dy,'m--') legend('原数据','Spline后','Spline求导')
效果:
堪称完美!
本回答被题主选为最佳回答 , 对您是否有帮助呢?解决 无用评论 打赏 举报 编辑记录
悬赏问题
- ¥20 机器学习能否像多层线性模型一样处理嵌套数据
- ¥20 西门子S7-Graph,S7-300,梯形图
- ¥50 用易语言http 访问不了网页
- ¥50 safari浏览器fetch提交数据后数据丢失问题
- ¥15 matlab不知道怎么改,求解答!!
- ¥15 永磁直线电机的电流环pi调不出来
- ¥15 用stata实现聚类的代码
- ¥15 请问paddlehub能支持移动端开发吗?在Android studio上该如何部署?
- ¥20 docker里部署springboot项目,访问不到扬声器
- ¥15 netty整合springboot之后自动重连失效