图形算法库研发（五）：非线性方程组求解

图形算法库中的非线性方程组求解需求

非线性方程组求解是图形算法库核心中的核心、曲线与曲线求交、曲线与曲面求交、曲面与曲面求交、空间距离检测、曲面求交内环检测等等等等。
总的来说，图形算法库对非线性方程组求解需求概况如下：

• 高效
  这就决定了基于当前技术，我们只能采用数值类算法，不能采用符号计算类的算法。
• 能求全部根
  丢根就意味着丢交点，就意味着后续处理（例如布尔运算、裁剪等等）拓扑结构会出问题，也就意味着最后结果会出错。
• 误差稳定性
  尤其是高次方的时候的稳定性，这对数值算法来说也是一个很大的挑战。
• 无穷多根能高效正确给出反馈
  重合问题是曲线或者曲面求交的重大难题，重合对于方程组来说就是无穷多组根，如何高效正确的给出反馈是个难点。
• 支持方程数和变量数不一致
  一般的解方程都要求一致，但某些情况下，例如曲线求交重合情况下，这两个个数不一致方便我们增加方程来求某些特征。

实现方案

基于上面的需求，WGP实现了一个方程组求解器，它的核心算法是在区间迭代法上做了改造。选择区间迭代法作为基础算法是因为它是目前已知的能稳定求出定义域内所有根的唯一有效数值算法。
当然，区间迭代法有很多不足之处，我们只是借鉴了它的思想，设计出了我们独有的数值算法。
关于区间迭代法思想的说明，可以看我下面的这篇文章。
数值法解非线性方程组全部实根最可行的算法——区间迭代法的几何意义

代码示例

WGP的方程组求解器，由一个方程组定义（EquationSystem）和一个求解器组成（Solver），下面是一个直接调用最基础源码的例子。

inline void TestSolver1() {
	class MyEquationSystem : public wgp::StandardEquationSystem {
	public:
		//get the count of equations
		virtual int GetEquationCount() {
			return 2;
		}
		//get the count of variables
		virtual int GetVariableCount() {
			return 2;
		}
		//get the tolerance of the ith variable
		virtual double GetVariableEpsilon(int i) {
			return 1E-7;
		}
		//get the tolerance of the ith equation
		virtual double GetValueEpsilon(int i) {
			return 1E-12;
		}
		//calculate values of equations
		virtual void CalculateValue(const wgp::StandardEquationsVariable& variable, wgp::StandardIntervalVector& value) {
			/*
			x^2 - y = 0
			x - y^2 = 0
			*/
			wgp::Interval x = variable.Get(0);
			wgp::Interval y = variable.Get(1);
			*value.Get(0) = sqr(x) - y;
			*value.Get(1) = x - sqr(y);
		}
		//calculate partial derivatives of equations
		virtual void CalculatePartialDerivative(const wgp::StandardEquationsVariable& variable, wgp::StandardIntervalMatrix& value) {
			/*
			2 * x,   -1
			1   ,    -2 * y
			*/
			wgp::Interval x = variable.Get(0);
			wgp::Interval y = variable.Get(1);
			*value.Get(0, 0) = 2 * x;
			*value.Get(0, 1) = -1;
			*value.Get(1, 0) = 1;
			*value.Get(1, 1) = -2 * y;
		}
	};
	MyEquationSystem* equation_system = new MyEquationSystem();
	wgp::Solver<MyEquationSystem, wgp::StandardEquationsVariable, wgp::StandardIntervalVector, wgp::StandardIntervalMatrix, wgp::StandardMatrix> solver;
	wgp::StandardEquationsVariable initial_value(2);
	initial_value.Set(0, wgp::Interval(-100, 100));	//x
	initial_value.Set(1, wgp::Interval(-100, 100));	//y
	solver.SetEquationSystem(equation_system);
	solver.SetInitialVariable(initial_value);
	const wgp::Array<wgp::StandardEquationsVariable>& clear_roots = solver.GetClearRoots();
	const wgp::Array<wgp::StandardEquationsVariable>& fuzzy_roots = solver.GetFuzzyRoots();
	delete equation_system;

	std::cout << "root count: " << clear_roots.GetCount() << "\n";
	for (int i = 0; i < clear_roots.GetCount(); ++i) {
		std::cout <<
			"x=" << clear_roots.GetPointer(i)->Get(0).Center() << ", " <<
			"y=" << clear_roots.GetPointer(i)->Get(1).Center() << "\n";
	}
	if (fuzzy_roots.GetCount() > 0) {
		std::cout << "fuzzy" << "\n";
	}
}

这段代码只是个最基础类的演示，显得有点繁琐，下面网站做了简单封装
http://www.fangchengzu.com/download
当然，封装的因为要编辑方程组和考虑通用性，效率肯定会比直接用原始的偏低。