在 GNU 科学图书馆多根查找器中选择起点
Choosing starting point in GNU Scientific Library multiroot finder
我正在使用 GNU 科学图书馆的多根查找器的实现来求解以下非线性方程组中的未知数(x
和 y
):
不过,我对“起点”有点困惑:
Solve(const double *x, int maxIter = 0, double absTol = 0, double relTol = 0)
Find the root starting from the point X; Use the
number of iteration and tolerance if given otherwise use default
parameter values which can be defined by the static method SetDefault
起点是如何选择的?
在一维多项式问题中,这是Real-root_isolation的well-studied域。在二维方面,它不太为人所知。
首先注意我们可以将方程式转化为多项式,取第一个
sqrt[(x-x1)^2 + (y-y1)^2] + s (t2 -t1) = sqrt[(x-x1)^2 + (y-y1)^2]
两边平方
[(x-x1)^2 + (y-y1)^2] + 2 sqrt[(x-x1)^2 + (y-y1)^2][s (t2 -t1)] + [s (t2 -t1)]^2 = [(x-x1)^2 + (y-y1)^2]
重新排列
2 sqrt[(x-x1)^2 + (y-y1)^2][s (t2 -t1)] = [(x-x1)^2 + (y-y1)^2] - [(x-x1)^2 + (y-y1)^2] - [s (t2 -t1)]^2
再平方
4 [(x-x1)^2 + (y-y1)^2] [s (t2 -t1)]^2 = ([(x-x1)^2 + (y-y1)^2] - [(x-x1)^2 + (y-y1)^2] - [s (t2 -t1)]^2)^2
给我们一个多项式。(注)
一旦我们有了一组多项式,您就可以使用一些代数技巧。
我使用的一种技术是将多项式转换为 Bernstein polynomials。首先,重新缩放您的域,使两个变量都位于 [0,1] 范围内。如果 m, n 是两个多项式的次数,那么一个多项式可以表示为 sum
sum i=0..m sum j=0..n b_ij mCi nCj x^i (1-x)^(n-i) y^j (1-y)^(n-j)
其中mCi、nCj为二项式系数
这些多项式有很好的凸性属性。如果系数 b_ij 全部为正,则多项式的值对于 [0,1] 中的所有 x,y 都将为正。如果系数均为负,则类似。
这让你可以说“这个区域没有解决方案”。
所以解决问题的策略可能是:
- 选择解可能出现的最大区域。
- 将其细分为一组更小的区域
- 对于每个区域,计算每个方程的伯恩斯坦多项式
- 检查伯恩斯坦多项式的系数,如果它们都是正的(或都是负的)则拒绝该区域
- 您现在应该已经拒绝了大部分域。使用每个剩余区域中的一个点启动多根查找器。
如果您想了解如何构造 Berstein 多项式的详细信息,可以在 A new method for drawing Algebraic Surfaces
阅读我的论文
注意:通过平方我们实际上得到了更多我们想要的解决方案。在最初的问题中,我们想要原则 sqareroot,即 +sqrt(A),也有使用另一个根 -sqrt(A) 的解决方案。这使得问题更容易一些,而不是我们从两个双曲线的交集得到的四个解决方案,我们只有一个分支的交集,所以问题的一个或两个解决方案。
对于您的问题,有一种非常简单的方法可以获取起点之一。
假设s=0。然后每个方程将给出与两点等距的点集。这是一条线,点连线的线段的垂直平分线,那么简单的求三个垂直平分线的交点。这将是 Circumscribed_circle 的中心。这对于算法来说可能就足够了。更简单的是简单地取三个点的平均值。
我正在使用 GNU 科学图书馆的多根查找器的实现来求解以下非线性方程组中的未知数(x
和 y
):
不过,我对“起点”有点困惑:
Solve(const double *x, int maxIter = 0, double absTol = 0, double relTol = 0)
Find the root starting from the point X; Use the number of iteration and tolerance if given otherwise use default parameter values which can be defined by the static method SetDefault
起点是如何选择的?
在一维多项式问题中,这是Real-root_isolation的well-studied域。在二维方面,它不太为人所知。
首先注意我们可以将方程式转化为多项式,取第一个
sqrt[(x-x1)^2 + (y-y1)^2] + s (t2 -t1) = sqrt[(x-x1)^2 + (y-y1)^2]
两边平方
[(x-x1)^2 + (y-y1)^2] + 2 sqrt[(x-x1)^2 + (y-y1)^2][s (t2 -t1)] + [s (t2 -t1)]^2 = [(x-x1)^2 + (y-y1)^2]
重新排列
2 sqrt[(x-x1)^2 + (y-y1)^2][s (t2 -t1)] = [(x-x1)^2 + (y-y1)^2] - [(x-x1)^2 + (y-y1)^2] - [s (t2 -t1)]^2
再平方
4 [(x-x1)^2 + (y-y1)^2] [s (t2 -t1)]^2 = ([(x-x1)^2 + (y-y1)^2] - [(x-x1)^2 + (y-y1)^2] - [s (t2 -t1)]^2)^2
给我们一个多项式。(注)
一旦我们有了一组多项式,您就可以使用一些代数技巧。
我使用的一种技术是将多项式转换为 Bernstein polynomials。首先,重新缩放您的域,使两个变量都位于 [0,1] 范围内。如果 m, n 是两个多项式的次数,那么一个多项式可以表示为 sum
sum i=0..m sum j=0..n b_ij mCi nCj x^i (1-x)^(n-i) y^j (1-y)^(n-j)
其中mCi、nCj为二项式系数
这些多项式有很好的凸性属性。如果系数 b_ij 全部为正,则多项式的值对于 [0,1] 中的所有 x,y 都将为正。如果系数均为负,则类似。
这让你可以说“这个区域没有解决方案”。
所以解决问题的策略可能是:
- 选择解可能出现的最大区域。
- 将其细分为一组更小的区域
- 对于每个区域,计算每个方程的伯恩斯坦多项式
- 检查伯恩斯坦多项式的系数,如果它们都是正的(或都是负的)则拒绝该区域
- 您现在应该已经拒绝了大部分域。使用每个剩余区域中的一个点启动多根查找器。
如果您想了解如何构造 Berstein 多项式的详细信息,可以在 A new method for drawing Algebraic Surfaces
阅读我的论文注意:通过平方我们实际上得到了更多我们想要的解决方案。在最初的问题中,我们想要原则 sqareroot,即 +sqrt(A),也有使用另一个根 -sqrt(A) 的解决方案。这使得问题更容易一些,而不是我们从两个双曲线的交集得到的四个解决方案,我们只有一个分支的交集,所以问题的一个或两个解决方案。
对于您的问题,有一种非常简单的方法可以获取起点之一。
假设s=0。然后每个方程将给出与两点等距的点集。这是一条线,点连线的线段的垂直平分线,那么简单的求三个垂直平分线的交点。这将是 Circumscribed_circle 的中心。这对于算法来说可能就足够了。更简单的是简单地取三个点的平均值。