清空筛选条件
清空筛选条件

curvature of a discrete function

49 次查看（过去 30 天）
David Kusnirak
David Kusnirak 2013-1-16
编辑： Jan 2022-3-13
Hello,
I need to compute a curvature of a simple 2D discrete function like this one:
x=1:0.5:20;
y=exp(x);
can anybody help how to do that? thanks

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

采纳的回答

Jan
Jan 2013-1-16
编辑：Jan 2022-3-13
Your function seems to be a 1D function.
Are you looking for the 2nd derivative? While diff calculates the one-sided differential quotient, gradient uses the two-sided inside the interval:
gradient(gradient(y))
If you mean the curvature as reciprocal radius of the local fitting circle:
dx = gradient(x);
ddx = gradient(dx);
dy = gradient(y);
ddy = gradient(dy);
num = dx .* ddy - ddx .* dy;
denom = dx .* dx + dy .* dy;
denom = sqrt(denom) .^ 3;
curvature = num ./ denom;
curvature(denom < 0) = NaN;
Please test this, because I'm not sure if I remember the formulas correctly.
  3 个评论
显示 1更早的评论
Jan
Jan 2013-1-16
编辑：Jan 2013-1-16
Therefore I'm using an efficient C-Mex function: FEX: DGradient, which is 10 to 20 times faster and handles unevenly spaced data more accurate.
Jan
Jan 2013-1-16
编辑：Jan 2022-3-13
x = rand(1, 1e6);
tic; ddx = gradient(gradient(x)); toc
tic; ddx = DGradient(DGradient(x)); toc
tic; ddx = conv(x,[.25 0 -.5 0 .25],'same'); toc
Elapsed time is 0.251547 seconds.
Elapsed time is 0.025728 seconds.
Elapsed time is 0.028208 seconds
Matlab 2009a/64, Core2Duo, Win7

请先登录，再进行评论。

更多回答（2 个）

Roger Stafford
Roger Stafford 2013-1-16
编辑：Bruno Luong 2022-3-13
Let (x1,y1), (x2,y2), and (x3,y3) be three successive points on your curve. The curvature of a circle drawn through them is simply four times the area of the triangle formed by the three points divided by the product of its three sides. Using the coordinates of the points this is given by:
K = 2*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)) ./ ...
sqrt(((x2-x1).^2+(y2-y1).^2).*((x3-x1).^2+(y3-y1).^2).*((x3-x2).^2+(y3-y2).^2));
You can consider this as an approximation to the curve's curvature at the middle point of the three points.
Matt J
Matt J 2013-1-16
编辑：Matt J 2013-1-16
diff(y,2)./0.5^2

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by