如何在Octave中判断点是否属于参数化花瓣曲线？

我有一个“花瓣”的参数方程：
$$
x(t) = (r+\cos(nlept))\cos(t) \
y(t) = (r+\cos(nlept))\sin(t)
$$
其中nlep为花瓣数量，r为半径。我想要生成一条从(0,0)出发、在该曲线内部反弹（反射）的轨迹。遇到的问题是无法判断一个点是否属于该曲线，而这是计算转向时机的必要条件。由于无法消去参数t，我尝试求解非线性方程组，但不确定如何设置求解时t的初始值。

以下是我的代码：

t = linspace(-10, 10, 100);
x = (r + cos(nlep*t)).*cos(t);
y = (r + cos(nlep*t)).*sin(t);
plot(x, y, 'g')
hold on

例如，我沿一条水平线行进，想要确定属于该曲线的点：

p = linspace(min(x), max(x), 100);
for j=1:length(p)
  f = @(t) ([((r + cos(nlep*t)).*cos(t))-p(j); ((r + cos(nlep*t)).*sin(t))-3]);
  [t0, fval, info] = fsolve(f, t(j));
  if info==1 % info is never 1
    plot(p(j), 3, 'r')  
  endif
endfor

但除曲线本身外没有任何内容被绘制出来。请问正确的实现方式是什么？能否提供一些提示？感谢。

1. 优化初始值，提高fsolve收敛概率 #

你的参数方程本质是极坐标形式：$\rho(t) = r + \cos(nlep \cdot t)$，$x = \rho(t)\cos(t)$，$y = \rho(t)\sin(t)$。当前用循环的t(j)作为初始值完全不匹配目标点，这是fsolve无法收敛的核心原因。

正确的初始值应基于目标点的极坐标推导：

• 计算目标点的极角：$\theta_0 = \arctan2(3, p_j)$
• 用$\theta_0$作为t的初始值，因为曲线在极角$\theta$附近的参数t接近$\theta$

修改后的代码片段：

for j=1:length(p)
  f = @(t) ([((r + cos(nlep*t)).*cos(t))-p(j); ((r + cos(nlep*t)).*sin(t))-3]);
  % 改用目标点的极角作为初始值
  t_init = atan2(3, p(j));
  [t0, fval, info] = fsolve(f, t_init);
  % 放宽收敛判断，只要残差足够小就认为找到解
  if info == 1 || norm(fval) < 1e-6
    plot(p(j), 3, 'r')  
  endif
endfor

2. 更高效的交点筛选方式 #

遍历横坐标找交点效率低且容易出错，直接遍历曲线的参数t，筛选出符合条件的点更可靠：

% 玫瑰线周期为2π，只需遍历0到2π即可
t = linspace(0, 2*pi, 1000);
x = (r + cos(nlep*t)).*cos(t);
y = (r + cos(nlep*t)).*sin(t);
plot(x, y, 'g')
hold on

% 筛选y坐标接近3的点（设置容差）
tol = 1e-3;
match_idx = abs(y - 3) < tol;
plot(x(match_idx), y(match_idx), 'r', 'Marker', 'o')

3. 反弹轨迹的核心实现思路 #

要生成内部反弹轨迹，核心是计算射线与曲线的交点，而非判断点是否在曲线上：

• 从当前点$(x_0,y_0)$沿方向向量$\vec{v}=(dx, dy)$发射射线，参数化为：$x = x_0 + k\cdot dx$，$y = y_0 + k\cdot dy$（$k>0$）
• 联立射线方程与花瓣曲线的参数方程，求解最小的$k>0$，得到交点$(x_i,y_i)$
• 计算曲线在交点处的单位法向量$\vec{n}$（指向曲线内部），用反射公式计算新方向：$\vec{v}{new} = \vec{v}{old} - 2(\vec{v}_{old} \cdot \vec{n})\vec{n}$
• 重复上述步骤，逐步生成轨迹

曲线的切向量与法向量计算：

• 切向量：$\vec{t} = (x'(t), y'(t)) = \left(-\rho(t)\sin(t) + \rho'(t)\cos(t), \rho(t)\cos(t) + \rho'(t)\sin(t)\right)$，其中$\rho'(t) = -nlep \cdot \sin(nlep \cdot t)$
• 法向量：将切向量逆时针旋转90度得到$\vec{n} = (-y'(t), x'(t))$，再判断其是否指向内部（可通过$\vec{n}$与$(x_i,y_i)$的点积符号判断，若点积为负则指向内部）

内容的提问来源于stack exchange，提问作者Snork Maiden