Octave中如何判断一个点是否属于参数曲线？

Question

我有一个“花”的参数方程：

x(t) = (r+cos(nlep*t))cos(t),

y(t) = (r+cos(nlep*t))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

但是除了曲线本身之外什么也没有绘制。应该用什么方法来做呢？有人可以给个提示吗？预先感谢。

Answer 1

您的系统通常没有解 - 它是一个由两个方程组成的系统，其中有一个未知数，

t

。

fsolve

只能最小化两个函数结果的范数。

如果我理解正确的话，你想找到曲线与直线的交点，比方说

y = a * x + b

，已知

a

和

b

。然后你必须解一个方程：

f = @(t) a * (r + cos(nlep*t)).*cos(t) + b - (r + cos(nlep*t)).*sin(t);

即

a * x + b - y

与

x

和

y

被曲线参数形式替代。

这是特定情况下的曲线、直线和解图：

r = 1; nlep = 5;

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

a = 1; b = 0.5; # y = x + 0.5

# plot the line
p = linspace(min(x), max(x), 100);
yp = a * p + b;
plot(p, yp, '-b')

# find and plot an intersection point
f = @(t) a * (r + cos(nlep*t)).*cos(t) + b - (r + cos(nlep*t)).*sin(t);
[t0, fval, info] = fsolve(f, t(j));

if info==1
  x0 = (r + cos(nlep*t0))*cos(t0);
  y0 = (r + cos(nlep*t0))*sin(t0);
  disp([x0, y0]);
  # plot the solution point
  plot(x0, y0, '*r', 'MarkerSize', 10)
endif

hold off

这里的问题是这个程序只能找到一个解。要找到许多/所有解决方案， 人们必须尝试使用不同的 t 初始值。这是最后部分的可能版本 上一段代码：

# find and plot all intersection points
x0 = [];
y0 = [];

f = @(t) a * (r + cos(nlep*t)).*cos(t) + b - (r + cos(nlep*t)).*sin(t);

for tinitial = linspace(-pi, pi, 20)
  [t0, fval, info] = fsolve(f, tinitial);
  if info==1
    x0new = (r + cos(nlep*t0))*cos(t0);
    y0new = (r + cos(nlep*t0))*sin(t0);
    #disp([x0i, y0i]);

    # check if the new solution wasn't already found
    exists = false;
    for i = 1:length(x0)
      if(abs(x0(i)-x0new) + abs(y0(i)-y0new) < 1e-5)
        exists = true;
        break;
      endif
    endfor

    if not(exists) #unefficient, but it happens only a few times.
      x0 = [x0, x0new];
      y0 = [y0, y0new];
    endif
  endif
endfor

disp(sprintf('%d solutions found', length(x0)))

# plot the solution points
plot(x0, y0, '*r', 'MarkerSize', 10)

hold off