natural cubic spline interpolation of y-values: how to get derivative of the spline wrt the y-values?

14 次查看（过去 30 天）
SA-W
SA-W 2023-3-8
编辑： Torsten 2023-3-9
Given a data set with support points x_1,...,x_n and corresponding y-values y_1,...,y_n.
My objective is to create a cubic spline f (with natural boundary conditions) that passes through the y_values. There are, of course, plenty of functions for doing this.
However, for a parameter identification procedure, I have to compute the derivative of the spline f with respect to the y-values -- at arbitrary points within [x1, x_n].
Is there an easy way using built-in functions of Matlab to compute the sensitivities?
  2 个评论
Torsten
Torsten 2023-3-8
I can understand that you have to compute the derivative of the spline with respect to the parameters, but why with respect to the y-values ?
SA-W
SA-W 2023-3-8
Because the values y_1,...,y_n are the parameters of my optimization. My objective is to optimize the values y_1,...,y_n and interpolate between them with cubic splines. Is that clear?

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2023-3-8
编辑：Bruno Luong 2023-3-8
Ran in:
The derivative f wrt to y_i is the spline interpolate b_i := (0,0,...,1,0...) where 1 is at ith position, since the spline is linear to y values.
Replace spline command with your function that computes natural spline pp form
x=cumsum(rand(1,10));
y=rand(size(x));
xi=linspace(min(x),max(x),500);
f=ppval(spline(x,y),xi)
f = 1×500
0.8578 0.8215 0.7861 0.7516 0.7180 0.6852 0.6533 0.6223 0.5920 0.5627 0.5341 0.5063 0.4794 0.4532 0.4278 0.4032 0.3794 0.3563 0.3340 0.3124 0.2915 0.2714 0.2519 0.2332 0.2152 0.1978 0.1812 0.1652 0.1498 0.1351
plot(x,y,'or',xi,f,'b') %
b=eye(length(x));
yd=spline(x,b);
dfdy=ppval(yd,xi); % dfdy(i,j) is the derivative of f(xi(j)) with respect to y(i))
figure
plot(xi,dfdy')
  17 个评论
显示 15更早的评论
SA-W
SA-W 2023-3-9
That said, if I wanted to calculate
dspline/dy = (spline(y + e_i * h) - spline(y)) / h
by finite differencing,
I could choose an arbitrary value for h as it cancels out anyway, right?

请先登录，再进行评论。

更多回答（2 个）

Bruno Luong
Bruno Luong 2023-3-8
编辑：Bruno Luong 2023-3-8
In this thread https://fr.mathworks.com/matlabcentral/answers/1894800-how-to-remove-noise-from-curves-and-take-their-derivates?s_tid=srchtitle
you can find my function that compute the derivative of a piecewise polynomiall function (pp), inclusing pp form of the spline functions. This function returns the pp form of the derivative, so you can evaluate using MATLAB ppval.
function ppd = ppder(pp)
ppd = pp;
coefs = ppd.coefs;
n = size(coefs,2);
ppd.coefs = coefs(:,1:n-1).*(n-1:-1:1);
ppd.order = ppd.order-1;
end
  4 个评论
显示 2更早的评论
SA-W
SA-W 2023-3-8
The derivative f wrt to y_i is the spline interpolate b_i := (0,0,...,1,0...) where 1 is at ith position, since the spline is linear to y values.
This is not clear to me. If I take the derivative of f (a scalar function) with respect to y_i (a scalar), the result should be a single number and not a vector b_i := (0,0,...,1,0...)...
Also, why is a cubic spline linear to y values? If I change yn to yn+h, f will be more affected at yn than at y1, for instance.
Bruno Luong
Bruno Luong 2023-3-8
编辑：Bruno Luong 2023-3-8
See my answer below, but
"the result should be a single number"
No the result is a scalar function. If we take f at a given point x then it is a scalar.
"...and not a vector b_i := (0,0,...,1,0...)."
I did not tell the derivative is b, the derivative is the spline interpolating b
"why is a cubic spline linear to y values?"
You clearly missunderstand and male confusion betwen being linear and being a lilnear function

请先登录，再进行评论。

Torsten
Torsten 2023-3-8
移动：Torsten 2023-3-8
Why do you want to compute the sensitivities manually ?
Usually, the fitting software computes them using a finite-difference approximation, i.e. by calling your function with
y_i
and
y_i+h
getting back
f_j(y_1,...,y_i,...,y_n) and f_j(y_1,...,y_i+h,...,y_n)
and approximating
df_j/dy_i = (f_j(y_1,...,y_i+h,...,y_n)-f_j(y_1,...,y_i,...,y_n))/h
And this would also be my suggestion on how to do it manually if it is really needed.

类别

在 Help Center 和 File Exchange 中查找有关 Splines 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by