四次方程根式解+四次以上方程近似解的js实现代码（下）——方程求解的核心类

跟范盛金大大不一样（我呸，我又厚颜无耻地跟大数学家比了），我在解一元高次方程方面没有加入自己的任何想法，也就是说，以下代码的实现全是抄现成的算法，因此不屑一提。如果一定要挖些亮点的话，也有，就是我会针对高次方程的一些特殊情况做了有针对性的降次，从而简化部分的运算以及减少精度的损失。如ax^4+cx^2+e=0这样的，我会先换元降到2次然后再求解。

废话不说了，类代码上完，本篇就结束了。感谢一直关注《矩阵的史诗级玩法》系列的粉丝们！

FormulaUtil.js

/**
 * 方程求解工具类
 * 由于二次或以上的方程涉及复数，因此所有的解都用复数存储
 * 
 */
function FormulaUtil()
{
	
}

/**
 * 计算方程ax+b=0的解 
 * 一次更没什么好说的了，小学的东西
 * @param a
 * @param b
 * @return 
 * 
 */
FormulaUtil.calcLineFormulaZero = function(a, b)
{
	var results = [];
	//最高次项系数为0，方程无解或者有无穷个解，此处忽略
	if(a != 0)
	{
		results.push(new ComplexNum(-b / a));
	}
	return results;
}

/**
 * 计算方程ax^2+bx+c=0的两个解 
 * 二次没什么好说的了，初中的东西
 * @param a
 * @param b
 * @param c
 * @return 
 * 
 */
FormulaUtil.calcSquareFormulaZero = function(a, b, c)
{
	//最高次项系数为0，用低一次的方程求解
	if(FormulaUtil.near0(a))
	{
		return FormulaUtil.calcLineFormulaZero(b, c);
	}
	var results = [];
	//二次判别式部分
	var deltaRoot = new ComplexNum(b * b - 4 * a * c).squareRoot();
	for(var i = 0; i < 2; i ++)
	{
		results.push(new ComplexNum(-b).add(deltaRoot[i]).divideNumber(2 * a));
		
	}
	return results;
}

/**
 * 计算方程ax^3+bx^2+cx+d=0的3个解 
 * 这里用一个我到现在都还没搞懂原理但形式足够简洁的盛金公式法来实现
 * @param a 3次项系数
 * @param b 2次项系数
 * @param c 1次项系数
 * @param d 常数项
 * @return 
 * 
 */
FormulaUtil.calcCubicFormulaZero = function(a, b, c, d)
{
	//最高次项系数为0，用低一次的方程求解
	if(FormulaUtil.near0(a))
	{
		return FormulaUtil.calcSquareFormulaZero(b, c, d);
	}
	
	//加入特殊情况处理，目的在于减少不必要的精度损失
	
	//形如ax^3+d=0的情况，直接开方即可
	if(FormulaUtil.near0(b) && FormulaUtil.near0(c))
	{
		return new ComplexNum(-d / a).cubicRoot();
	}
	
	//ax^3+bx^2=0，常数项和一次项都为0，两个根为0，然后降为一次即可
	if(FormulaUtil.near0(d) && FormulaUtil.near0(c))
	{
		var lineRoots = FormulaUtil.calcLineFormulaZero(a, b);
		lineRoots.unshift(new ComplexNum());
		lineRoots.unshift(new ComplexNum());
		return lineRoots;
	}
	
	//ax^3+bx^2+cx=0，常数项为0，一个根为0，然后降为两次即可
	if(FormulaUtil.near0(d))
	{
		var squareRoots = FormulaUtil.calcSquareFormulaZero(a, b, c);
		squareRoots.unshift(new ComplexNum());
		return squareRoots;
	}			
	
	var A = b * b - 3 * a * c;
	var B = b * c - 9 * a * d;
	var C = c * c - 3 * b * d;
	
	var DELTA = B * B - 4 * A * C;
	
	var i = 0;
	var results = [];
	
	var x1;
	var x;
	
	console.log(DELTA);
	
	if(A == 0 && B == 0)
	{
		for(i = 0; i < 3; i ++)
		{
			results.push(new ComplexNum(-b / 3 / a));
		}
	}else if(DELTA > 0)
	{
		//这里本来想偷懒下，用二次方程的公式套进去，但没想到A=0就不能用了
		var roots = [];
		//二次判别式部分
		var deltaRoot = new ComplexNum(B * B - 4 * A * C).squareRoot();
		for(i = 0; i < 2; i ++)
		{
			roots.push(new ComplexNum(-B).add(deltaRoot[i]).divideNumber(2));
		}
		var ys = [];
		//y的立方根
		var ycrs = [];
		for(i = 0, len = roots.length; i < len; i ++)
		{
			var root = roots[i];
			var y = root.multiplyNumber(a * 3).addNumber(A * b);
			ys.push(y);
			ycrs.push(new ComplexNum((y.real > 0 ? 1 : -1) * Math.pow(Math.abs(y.real), 1 / 3)));
		}
		x1 = new ComplexNum(-b).subtract(ycrs[0]).subtract(ycrs[1]).divideNumber(3 * a);
		console.log(ycrs[0], ycrs[1]);
		results.push(x1);
		for(i = 0; i < 2; i ++)
		{
			x = new ComplexNum(-2 * b).add(ycrs[0]).add(ycrs[1]).
				add(ycrs[0].subtract(ycrs[1]).multiply(ComplexConsts.I).
				multiplyNumber(i == 0 ? Math.sqrt(3) : -Math.sqrt(3))).divideNumber(6 * a);
			results.push(x);
		}				
	}else if(DELTA == 0)
	{
		var K = B / A;
		results.push(new ComplexNum(- b / a + K), - 1 / 2 * K, - 1 / 2 * K);				
	}else if(DELTA < 0)
	{
		//盛金公式上说这种情况下A必大于0，所以就没用复数来搞了
		var T = (2 * A * b - 3 * a * B) / 2 * Math.sqrt(A * A * A);
		var theta = Math.acos(T);
		x1 = new ComplexNum((-b - 2 * Math.sqrt(A) * Math.cos(theta / 3)) / 3 / a);
		for(i = 0; i < 2; i ++)
		{
			x = new ComplexNum(-b + Math.sqrt(A) * (Math.cos(theta / 3) + (i == 0 ? 1 : -1) * Math.sin(theta / 3) * Math.sqrt(3)) / 3 / a);
		}
	}			
	return results;
}

/**
 * 计算方程ax^4+bx^3+cx^2+dx+e=0的4个解 
 * 以下代码按照费拉里公式进行实现
 * @param a 4次项系数
 * @param b 3次项系数
 * @param c 2次项系数
 * @param d 1次项系数
 * @param e 常数项
 * @return 
 * 
 */
FormulaUtil.calcFourFormulaZero = function(a, b, c, d, e)
{
	
	//最高次项系数为0，用低一次的方程求解
	if(FormulaUtil.near0(a))
	{
		return FormulaUtil.calcCubicFormulaZero(b, c, d, e);
	}
	
	//测试发现，费拉里法无法处理一些特殊情况，m=0时出错几率较高，但未能找出必现规律，为此，对于一些能用简单方法直接降次数的情况，此处加入特殊处理
	//此外，特殊情况处理的好处是运算量小，可以减少中途计算中的精度损失
	
	//形如ax^4+e=0的情况，直接开方即可
	if(FormulaUtil.near0(b) && FormulaUtil.near0(c) && FormulaUtil.near0(d))
	{				
		return new ComplexNum(- e / a).getRoot(4);
	}		
	
	//ax^4+bx^3=0，常数项，一次项，二次项均为0，则有三个根为0，然后直接降到1次即可
	if(FormulaUtil.near0(c) == 0 && FormulaUtil.near0(d) && FormulaUtil.near0(e))
	{
		var lineRoots = calcLineFormulaZero(a, b);
		lineRoots.unshift(new ComplexNum());
		lineRoots.unshift(new ComplexNum());
		lineRoots.unshift(new ComplexNum());
		return lineRoots;
	}
	
	//ax^4+bx^3+cx^2=0，常数项和一次项为0，则有两个根为0，然后直接降到2次即可
	if(FormulaUtil.near0(d) && FormulaUtil.near0(e))
	{
		var squareRoots = FormulaUtil.calcSquareFormulaZero(a, b, c);
		squareRoots.unshift(new ComplexNum());
		squareRoots.unshift(new ComplexNum());
		return squareRoots;
	}		
	
	//ax^4+bx^3+cx^2+dx=0，常数项为0，有一个根为0，然后直接降到3次即可
	if(FormulaUtil.near0(e))
	{
		var cubicRoots = FormulaUtil.calcCubicFormulaZero(a, b, c, d);
		cubicRoots.unshift(new ComplexNum());
		return cubicRoots;
	}
	
	//ax^4+cx^2+e=0，一次项和三次项为0，为双二次方程，用换元法让y=x^2即可降为两次
	if(FormulaUtil.near0(b) && FormulaUtil.near0(d))
	{
		var twoSquareRoots = FormulaUtil.calcSquareFormulaZero(a, c, e);
		var fourRoots = twoSquareRoots[0].squareRoot().concat(twoSquareRoots[1].squareRoot());
		return fourRoots;
	}
	
	var P = (c * c + 12 * a * e - 3 * b * d) / 9;
	var Q = (27 * a * d * d + 2 * c * c * c + 27 * b * b * e - 72 * a * c * e - 9 * b * c * d) / 54;
	var D = new ComplexNum(Q * Q - P * P * P).squareRoot()[0];
	
	var u1 = new ComplexNum(Q).add(D).cubicRoot()[0];
	var u2 = new ComplexNum(Q).subtract(D).cubicRoot()[0];
	
	var u = u1.length > u2.length ? u1 : u2;			
	
	var v = (u.real == 0 && u.image == 0) ? new ComplexNum() :(new ComplexNum(P).divide(u));			
	
	var omegas = [ComplexConsts.OMEGA_ZERO, ComplexConsts.OMEGA, ComplexConsts.OMEGA_SQUARE, ComplexConsts.OMEGA_CUBIC];
	
	var currentMLength = 0;
	var maxM;
	var currentOmegaObj;
	//算出所有的m来，并取出最大值
	for(var k = 0; k < 3; k ++)
	{
		var m = new ComplexNum(b * b - 8 / 3 * a * c);
		var omegaObj = u.clone();//omegas[k].multiply(u);
		omegaObj = omegaObj.add(v);
		omegaObj = omegaObj.multiplyNumber(4 * a);
		m = m.add(omegaObj);	
		m = m.squareRoot()[0];
		if(m.length >= currentMLength)
		{
			currentOmegaObj = omegaObj;
			currentMLength = m.length;
			maxM = m;
		}
	}
	
	m = maxM;
	omegaObj = currentOmegaObj;
	
	var S = new ComplexNum(2 * b * b - 16 / 3 * a * c);
	S = S.subtract(omegaObj);
	
	var T;
	if(!FormulaUtil.near0(m.real) || !FormulaUtil.near0(m.image))
	{
		T = new ComplexNum(8 * a * b * c - 16 * a * a * d - 2 * b * b * b); 
		T = T.divide(m);
	}else
	{
		T = new ComplexNum();
	}
	
	var results = [];
	
	for(var n = 1; n <= 4; n ++)
	{
		var value = new ComplexNum();
		//-1的[n/2]次方，[n/2]表示小于n/2的最大整数
		var minus1Pow0_5 = n > 2 ? -1 : 1;
		value = value.subtractNumber(b);		
		value = value.add((new ComplexNum(minus1Pow0_5, 0).multiply(m)));
		value = value.add(S.add(new ComplexNum(minus1Pow0_5, 0).multiply(T)).squareRoot()[0].multiplyNumber(n % 2 == 0 ? -1 : 1));
		value = value.divideNumber(4 * a);
		results.push(value);
	}	
	
	return results;
}

/**
 * 计算多项式方程的解 
 * @param args
 * 
 */
FormulaUtil.calcPolyFormulaZero = function()
{		
	
	var currentArgs = arguments.concat();
	
	while(currentArgs.length > 0 && currentArgs[0] == 0)
	{
		currentArgs.shift();
	}
	
	var argCount = currentArgs.length;
	if(argCount == 0)
	{
		throw new Error("必须传入参数");
		return null;
	}
	
	if(argCount <= 1)
	{
		return [];
	}
	
	if(argCount <= 5)
	{
		return _polyZeroFunctions_vec[argCount - 1].apply(null, currentArgs);
	}
	var i, len;
	
	
	//导数后的系数
	//对于5次或以上的情况，可以通过导数降低次数，然后用切线法求出近似值
	var dArgs = FormulaUtil.getDArgs(currentArgs);
				
	//二阶导数的系数
	var d2Args = FormulaUtil.getDArgs(dArgs);
	
	//先算出导数等于0的点（可能是极值点）
	var dSolutions = FormulaUtil.calcPolyFormulaZero.apply(null, dArgs);	
	//解方程的过程可能会出现虚根，求近似值只能去掉了（暂没想到比较好的求近似虚根的方法）
	var dRealSolutions = FormulaUtil.getRealSolutions(dSolutions);
	dRealSolutions.sort(Array.NUMERIC);
	var yValues = [];
	var realSolutions = [];
	for(i = 0, len = dRealSolutions.length; i < len; i ++)
	{				
		var solution = dRealSolutions[i];
		yValues.push(getPolyValue(solution, currentArgs));
	}
	//无解的话可能是单调递增或者递减，这时候应该补上一个解，另外，最左侧的零点左侧和最右侧的零点右侧都可能有解，这些逻辑暂时还没做
	for(i = 0, len = yValues.length - 1; i < len; i ++)
	{
		if(yValues[i] * yValues[i + 1] < 0)
		{
			var xValue = findRoot(dRealSolutions[i], dRealSolutions[i + 1], currentArgs, dArgs, d2Args);
			console.log("在" + dRealSolutions[i] + "和" + dRealSolutions[i + 1] + "之间有一个解x=" + xValue);
			realSolutions.push(new ComplexNum(xValue));
		}else if(yValues[i] == 0)
		{
			realSolutions.push(new ComplexNum(dRealSolutions[i]));
		}else if(yValues[i + 1] == 0)
		{
			realSolutions.push(new ComplexNum(dRealSolutions[i + 1]));
		}
	}
	
	return realSolutions;
}

FormulaUtil.getPolyValue = function(xValue, args)
{
	var argCount = args.length;
	var value = 0;
	//遍历所有次方的项，然后累加
	for(var i = 0, len = argCount; i < len; i ++)
	{
		var currentK = args[i];
		//求出当前项的次数
		var times = argCount - 1 - i;
		value += currentK * (times == 0 ? 1 : Math.pow(xValue, times));
	}
	return value;
}

/**
 * 在指定区间内找到一个高次方程的解 
 * @param min 区间的下限
 * @param max 区间的上限
 * @param formulaArgs 高次方程的系数
 * @param dArgs 方程一阶导数的系数
 * @param d2Args 方程二阶导数的系数
 * @param precision 允许的最大误差
 * @return 
 * 
 */
FormulaUtil.findRoot = function(min, max, formulaArgs, dArgs, d2Args, precision)
{
	if(isNaN(precision))
	{
		precision = 0.00001;
	}
	//最好用牛顿切线法，但现在先用二分法凑合一下吧			
	do
	{
		var center = (min + max) / 2;			
		var minYValue = getPolyValue(min, formulaArgs);
		var maxYValue = getPolyValue(max, formulaArgs);
		var centerYValue = getPolyValue(center, formulaArgs);		
		if(centerYValue == 0)
		{
			return center;
		}
		if(centerYValue * minYValue < 0)
		{
			max = center;
		}else
		{
			min = center;
		}
	}while(Math.abs(min - max) >= precision);
	return center;
}

FormulaUtil.getDArgs = function(sourceArgs)
{
	var dArgs = [];
	var argCount = sourceArgs.length;
	for(var i = 0, len = argCount - 1; i < len; i ++)
	{
		var newK = (argCount - 1 - i);
		var finalK = newK * sourceArgs[i];
		dArgs.push(finalK);
	}
	return dArgs;
}

/**
 * 为了实用上的方便，这里提供一个通过复数解得到实数解的方法，可能有相等的实根 
 * @param complexNums
 * @param tolerance
 * @return 
 * 
 */
FormulaUtil.getRealSolutions = function(complexNums, tolerance)
{
	if(isNaN(tolerance))
	{
		tolerance = 0.00001;
	}
	var numbers = [];
	for(var i = 0, len = complexNums.length; i < len; i ++)
	{
		var item = complexNums[i];
		if(FormulaUtil.near0(item.image, tolerance))
		{
			numbers.push(item.real);
		}
	}
	return numbers;
}

FormulaUtil.near0 = function(num, tolerance)
{
	if(isNaN(tolerance))
	{
		tolerance = 0.0000001;
	}
	return Math.abs(num) <= tolerance;
}

//不同次数方程的求解方法列表
FormulaUtil._polyZeroFunctions = [FormulaUtil.calcLineFormulaZero, FormulaUtil.calcLineFormulaZero, FormulaUtil.calcSquareFormulaZero, FormulaUtil.calcCubicFormulaZero, FormulaUtil.calcFourFormulaZero];

白白了，接下来我会回到矩阵史诗级玩法的连载中，敬请继续关注！