请教一个关于BSplineFunction的问题
mathematica吧
全部回复
仅看楼主
level 3
駩爆🔥 楼主
接触mathematica没多久,最近想用B样条函数得到我需要的曲线方程,但BSplineFunction生成的函数是横纵坐标值关于t的函数,(t取0到1),我想得到图中y关于x的函数,有什么办法么?
另外,想求图上这个曲线在不同位置的曲率,应该怎么写?
pts1={{0,100},{200,100},{200,0}};
fend=BSplineFunction[pts1];
Show[Graphics[{Red,PointSize[0.03],Point[pts1],Green,Line[pts1]},Axes->True,AxesOrigin->{0,0},AxesLabel->{"x(mm)","y(mm)"},AxesStyle->Directive[Black,10]],ParametricPlot[fend[t],{t,0,1}]]
Plot[N[ArcCurvature[{fend[[1]][t],fend[[2]][t]},t]/.t->x],{x,0,1}]
2023年02月23日 00点02分 1
吧务
level 10
确定自变量范围和函数单值性的情况下,用Rescale写一个吧,改样条函数对象的表达式也行。
2023年02月23日 14点02分 2
吧务
level 10
fendx[t_?NumericQ] := fend[t][[1]] /; 0<=t<=1
fendy[t_?NumericQ] := fend[t][[2]] /; 0<=t<=1
Plot[ArcCurvature[{fendx@t,fendy@t},t]/.t->x,{x,0,1}]
2023年02月23日 14点02分 3
不行的,要算曲率就不能用黑箱函数,否则会触发Derivative的一个坑爹特性,参SE帖子《Why does Arg'[1. + I] return -0.5?》(编号196998)
2023年03月04日 05点03分
另一个相关帖:《Symbolic derivatives are being calculated numerically》(编号29329)
2023年03月04日 07点03分
在9楼又补了几句。
2023年03月04日 10点03分
level 3
駩爆🔥 楼主
@asdasd1dsadsa
非常感谢您的回答,第二个问题(也就是曲率的问题)应该已经解决了。
但第一个问题,用Rescale似乎还有些问题。您的意思是不是这样:
fendyx[x_] := fend[Rescale[x, {0, 200}]][[2]];
Plot[fendyx[x], {x, 0, 200}]
把x先转化为0到1的范围,然后再代入到b样条函数里求y,但问题是b样条函数的自变量t取0到1,和x取0到200,它们之间不是线性变化的关系。。即t取0.5的时候,x并不是100。
所以好像得出的曲线和b样条函数的曲线不一样。
我这个问题确定是自变量范围和函数,我就是用来解决一个绕制线圈的具体工程问题的。
2023年02月24日 11点02分 4
噢,确实。那可以试试自己对fend重新采样插值,做一个xy关系函数。
2023年02月24日 13点02分
level 3
駩爆🔥 楼主
@asdasd1dsadsa
谢谢,我按照您说的,在曲线上找50个点{x,y},然后进行样条插值,的确得到了和原B样条曲线看起来一样的曲线。
pts1 = {{0, 100}, {200, 100}, {200, 0}};
fend = BSplineFunction[pts1];
Show[Graphics[{Red, PointSize[0.03], Point[pts1], Green, Line[pts1]},
Axes -> True, AxesOrigin -> {0, 0}, AxesLabel -> {"x(mm)", "y(mm)"},
AxesStyle -> Directive[Black, 10]],
ParametricPlot[fend[t], {t, 0, 1}]]
numberPoints = 50;
interPoints =
Table[{fend[(i - 1)/numberPoints][[1]],
fend[(i - 1)/numberPoints][[2]]}, {i, numberPoints + 1}];
fendyx = Interpolation[ interPoints, Method -> "Spline"];
Show[Graphics[{Red, PointSize[0.03], Point[pts1], Green, Line[pts1]},
Axes -> True, AxesOrigin -> {0, 0}, AxesLabel -> {"x(mm)", "y(mm)"},
AxesStyle -> Directive[Black, 10]],
ParametricPlot[fend[t], {t, 0, 1}],
Plot[fendyx[x], {x, 0, 200},
PlotStyle -> Directive[Red, Dashing[Small]]]]
但求曲率又出现了新问题:
我又看了一下您3楼的求曲率的代码,画出曲率和t的图发现,t的取值范围并不是0到1,而仅仅有0.4至0.6附近,我加上PlotRange -> Full也没有变化,不知这是为什么。
另外,我还是希望得到曲率跟x的关系曲线,我用插值得到的函数求曲率,但得到的曲率在x接近200附近出现了振动,感觉很奇怪,B样条本身在这个地方的曲率看着不像能振荡的,不知道是不是插值函数的原因。。
2023年02月25日 08点02分 6
最后一幅图会振荡主要还是numberPoints小了。你取个1000试试。不过右端点的性质可能还要再仔细分析下。
2023年03月04日 05点03分
右侧是个可去间断点,参8楼。
2023年03月04日 05点03分
3楼代码算曲率时为什么会有问题我已在该楼楼中楼做了说明。
2023年03月04日 06点03分
吧务
level 15

用BSplineBasis可能更优一些。注意BSplineBasis可以用PiecewiseExpand展开:
《Convert BSplineFunction into two Interpolating Functions》
https://mathematica.stackexchange.com/q/19229/1871
2023年03月04日 04点03分 7
吧务
level 15
试着做了一下。对于这个只有3个点的问题,解析解是可解的。进一步参考
《What do the arguments of a generated BSplineFunction mean?》
https://mathematica.stackexchange.com/a/280442/1871
里的讨论,可以写出:
pts = {{0, 100}, {200, 100}, {200, 0}};
bs = BSplineFunction[pts];
m = Length[pts];
(* J.M. 的那个定义degree和knots的方法在这里不能用,因为点太少了,BSplineFunction会自动降阶 *)
degree = bs["Degree"][[1]];
knots = Rationalize@bs["Knots"][[1]];
{xu, yu} = Transpose[pts];
f[t_] = xu . Table[BSplineBasis[{degree, knots}, i - 1, t], {i, m}];
g[t_] = yu . Table[BSplineBasis[{degree, knots}, i - 1, t], {i, m}];
bs2[t_] = Simplify[{f[t], g[t]} // Map@PiecewiseExpand, 0 < t < 1]
func[x_] =
bs2[t][[2]] /. First@Solve[{bs2[t][[1]] == x, 0 < t < 1, x ∈ Reals}, t]
(* ConditionalExpression[-100 (-1 + (1 - Sqrt[200 - x]/(10 Sqrt[2]))^2),
0 < x < 200] *)
arc = ArcCurvature[{x, func[x]}, x] // FullSimplify
(* ConditionalExpression[(4 Sqrt[2/5])/(240 - 4 Sqrt[400 - 2 x] - x)^(3/2),
0 < x < 200] *)
(*对于0 < x < 200,x == 200处是个可去间断点*)
Limit[arc, x -> 200, Direction -> "FromBelow"]
(* 1/100 *)
2023年03月04日 05点03分 8
“可去间断点”这词好像用错了?应该叫“该点处的左极限存在”?
2023年03月04日 07点03分
吧务
level 15
3楼的楼中楼给的帖子可能还是不够直接,这里补个更易懂的例子:
Clear@func
func[x_?NumericQ] :=
Piecewise[{{Sin[x], 0 <= x <= 1}}, Echo["这个数被喂进来了鸭:" <> ToString@x]]
func'[0.4]
也就是说,Derivative的这个半吊子数值求导机制(使用高阶的中心差分公式对导数做数值计算)在撞到非数值量的时候会选择返回输入。这个其实倒也不难规避,真正的问题是用这个方法算出来的导数精度非常堪忧,所以,一定要避免使用Derivative算黑箱函数的数值导数。
2023年03月04日 10点03分 9
啊补充一句,InterpolatingFunction对符号求导有特殊的处理机制,所以,直接对InterpolatingFunction使用Derivative不会出现上述的问题。(这样会算出一个新的InterpolatingFunction。)
2023年03月04日 11点03分
BSplineFunction也是如此。
2023年03月04日 11点03分
那怎么搞?手动差分?
2023年03月04日 13点03分
@asdasd1dsadsa 这还真是个办法。NDSolve`FiniteDifferenceDerivative还挺好使的,某些场合这甚至是唯一可行的办法,参《Strange chattering in derivative of InterpolatingFunction returned by NDSolve using FEM》(SE编号276180)
2023年03月04日 13点03分
level 3
駩爆🔥 楼主
非常感谢吧主如此详细的解答,受益匪浅,果然6楼那个numberPoints改成1000曲率不会振荡了,算是解决了我实际的工程设计需求。
看了您其他楼层的回复,作为一一名初学者,也学到很多,感谢感谢
前几天有点忙,忘了及时回复了,抱歉。
2023年03月09日 09点03分 10
1