1. 程式人生 > >四次方程根式解+四次以上方程近似解的js實現程式碼(1)——複數類+複數常量+三角函式簡表

四次方程根式解+四次以上方程近似解的js實現程式碼(1)——複數類+複數常量+三角函式簡表

本人正在寫矩陣史詩級玩法系列部落格,寫到求二元二次方程組的地方來了,消元后最高會生成一元四次方程,而這個求根公式雖然成熟,但程式碼量也不少,所以單獨封裝成工具類。

本不打算講解的,但考慮到有的朋友可能沒接觸過複數,或者說雖然接觸過複數但已經忘得一乾二淨,那這裡我就簡單說一下好了。在實數範圍內,負數是沒辦法開平方的,而且對於一些簡單的曲線求交(二次或以下)問題來說,出現要對負數開平方的話就意味著沒有交點。然而在一些複雜的方程求解中,負數開平方也許是一個重要的中間過程,走完中間過程之後可能會得到實數根,所以數學家們定義了一個新的數i,名為虛數單位,並規定i^2=-1,同時可以使其參與到跟實數的四則運算中且交換律結合律等仍然成立。然後,實數1+實數2*i就構成了複數a+bi。

四次方程求解經常用到複數,其結構為實數部分一個係數(實部),虛數部分一個係數(虛部),程式碼如下。

function ComplexNum(real, image)
{
	this.real = isNaN(real) ? 0 : real;
	this.image = isNaN(image) ? 0 : image;
}

複數本身有四則運算,所以完整的程式碼我把這些運算的實現都寫進去了,加減最好理解,乘法稍有難度,除法和開方需要套路,如果覺得看不懂,那可以自己百度或者查閱相關的高中數學教材來補補這方面的知識。

ComplexNum.js

function ComplexNum(real, image)
{
	this.real = isNaN(real) ? 0 : real;
	this.image = isNaN(image) ? 0 : image;
	
	/**
	 * 設定複數的值 
	 * @param reala 複數的實部
	 * @param imagea 複數的虛部
	 * 
	 */
	function setTo(reala, imagea)
	{
		this.real = reala;
		this.image = imagea;
	}
	this.setTo = setTo;
	
	/**
	 * 複數相加 
	 * @param complex 要與其相加的複數
	 * @return 
	 * 
	 */
	function add(complex)
	{
		return new ComplexNum(this.real + complex.real, this.image + complex.image);
	}
	this.add = add;
	
	/**
	 * 為了使用上方便一點,此處加上與實數相加的方法 
	 * @param num 要與其相加的數字
	 * @return 
	 * 
	 */
	function addNumber(num)
	{
		return new ComplexNum(this.real + num, this.image + num);
	}
	this.addNumber = addNumber;
	
	/**
	 * 複數相減
	 * @param complex 要與其相減的複數
	 * @return 
	 * 
	 */
	function subtract(complex)
	{
		return new ComplexNum(this.real - complex.real, this.image - complex.image);
	}
	this.subtract = subtract;
	
	/**
	 * 為了使用上方便一點,此處加上與實數相減的方法 
	 * @param num 要與其相減的數字
	 * @return 
	 * 
	 */
	function subtractNumber(num)
	{
		return new ComplexNum(this.real - num, image);
	}
	this.subtractNumber = subtractNumber;
	
	/**
	 * 複數乘法 (運演算法則是利用多項式乘法進行展開得到)
	 * @param complexNum 要與其相乘的複數
	 * @return 
	 * 
	 */
	function multiply(complexNum)
	{
		return new ComplexNum(this.real * complexNum.real - this.image * complexNum.image, this.image * complexNum.real + this.real * complexNum.image);
	}
	this.multiply = multiply;
	
	/**
	 * 為了使用上方便一點,此處加上與實數相乘的方法 
	 * @param num 要與其相乘的數字
	 * @return 
	 * 
	 */
	function multiplyNumber(num)
	{
		return new ComplexNum(this.real * num, this.image * num);
	}
	this.multiplyNumber = multiplyNumber;
	
	/**
	 * 複數除法 (運演算法則是通過平方差公式,分子分母同時乘以分母的共軛複數以去除分母中的虛部,然後就利用乘法法則進行計算)
	 * @param complexNum 要與其相除的複數
	 * @return 
	 * 
	 */
	function divide(n)
	{
		//分母化為實數
		var denominator = n.real * n.real + n.image * n.image;
		//分子也乘以同樣的複數,併除以分母即得最終結果
		return new ComplexNum((this.real * n.real + this.image * n.image) / denominator, (this.image * n.real - this.real * n.image) / denominator);
	}
	this.divide = divide;
	
	/**
	 * 為了使用上方便一點,此處加上與實數相除的方法 
	 * @param num 要與其相除的數字
	 * @return 
	 * 
	 */
	function divideNumber(num)
	{
		return new ComplexNum(this.real / num, this.image / num);
	}
	this.divideNumber = divideNumber;
	
	/**
	 * 開平方運算
	 * @return 
	 * 
	 */
	function squareRoot()
	{			
		return this.getRoot(2);
	}
	this.squareRoot = squareRoot;
	
	/**
	 * 開立方運算 
	 * @return 
	 * 
	 */
	function cubicRoot()
	{
		return this.getRoot(3);
	}
	this.cubicRoot = cubicRoot;
	
	/**
	 * 開任意整數次方的運算 
	 * @param times
	 * @return 
	 * 
	 */
	function getRoot(times)
	{
		var vec = [];			
		//複數開方運算的原理是把輻角根據次數進行平分
		var degree = this.degree;
		degree /= times;
		//然後多個方根平分360度,所以需要算出每個方根之間的輻角間隔
		var degreeUnit = 360 / times;
		//複數長度(模)直接開方即可
		var lengthRoot = Math.pow(this.length, 1 / times);
		var cosDic = AngleUtil.getCosDic();
		var sinDic = AngleUtil.getSinDic();
		//然後就能通過迴圈生成所有開方結果
		for(var i = 0; i < times; i ++)
		{
			var currentDegree = (degree + i * degreeUnit);
			var currentAngle = currentDegree * Math.PI / 180;
			var cos = isNaN(cosDic[currentDegree]) ? Math.cos(currentAngle) : cosDic[currentDegree];
			var sin = isNaN(sinDic[currentDegree]) ? Math.sin(currentAngle) : sinDic[currentDegree];
			//trace(lengthRoot * cos, lengthRoot * sin);
			vec.push(new ComplexNum(lengthRoot * cos, lengthRoot * sin));
		}
		return vec;
	}
	this.getRoot = getRoot;
	
	/**
	 * 複數的輻角主值(複數所在點與座標連線跟水平線的夾角),以弧度為單位
	 * 必要時可用degree屬性可以更有效避免邊緣位置的浮點誤差引發的錯誤
	 * @see degree()
	 * 
	 */	
	Object.defineProperty(this, "angle", {get: getAngle});
	function getAngle()
	{
		//由於複數開方基於角度,所以一旦有浮點誤差就會導致角度在邊緣位置出現突躍導致嚴重錯誤,所以遇到浮點誤差的話,我會強制將角度設定為0
		if(Math.abs(this.image) < 0.0000001)
		{
			return this.real > 0 ? 0 : Math.PI; 
		}else if(Math.abs(this.real) < 0.0000001)
		{
			return this.image > 0 ? Math.PI * 0.5 : (this.image == 0 ? 0 : - Math.PI * 0.5);
		}else
		{
			return Math.atan2(this.image, this.real);
		}
	}

	
	/**
	 * 複數的輻角主值,以角度值表示(跟弧度相比,該方法能更有效地避免邊緣位置的浮點誤差引發的錯誤)
	 * 
	 */	
	Object.defineProperty(this, "degree", {get: getDegree});	 
	function getDegree()
	{
		//由於複數開方基於角度,所以一旦有浮點誤差就會導致角度在邊緣位置出現突躍導致嚴重錯誤,所以遇到浮點誤差的話,我會強制將角度設定為0
		if(Math.abs(this.image) < 0.0000001)
		{
			return this.real > 0 ? 0 : 180; 
		}else if(Math.abs(this.real) < 0.0000001)
		{
			return this.image > 0 ? 90 : (this.image == 0 ? 0 : 270);
		}else
		{
			return Math.atan2(this.image, this.real) / Math.PI * 180;
		}
	}
	
	/**
	 * 複數的模(長度) 
	 * @return 
	 * 
	 */
	Object.defineProperty(this, "length", {get: getLength});
	function getLength()
	{
		return Math.sqrt(this.real * this.real + this.image * this.image);
	}
	
	/**
	 * 返回當前複數的一個副本 
	 * @return 
	 * 
	 */
	function clone()
	{
		return new ComplexNum(this.real, this.image);
	}
	this.clone = clone;
	
	function toString()
	{
		var realStr = String((this.real != 0 || this.image == 0) ? real : "");
		var imageStr = (this.image == 0) ? "" : ((this.image < 0 ? ("-" + (-this.image)) : ((realStr == "" ? "" : "+") + this.image)) + "i");
		return realStr + imageStr;
	}
	this.toString = toString;
}

複數乘法跟三角函式密切相關,而且60度倍數的三角函數出現頻率特別高,所以我又定義了一些常量,以防直接計算累積不必要的浮點誤差。

ComplexConsts.js

function ComplexConsts()
{
	
}

/**
 * OMEGA是模為1,幅角主值等於120度的複數,它是1開三次方的結果,在解3次和4次方程中非常常用 
 * OMEGA就是希臘字母裡最像w的那個
 * 
 */
ComplexConsts.OMEGA = new ComplexNum(1 / 2, Math.sqrt(3) / 2);

/**
 * 理論上OMEGA的平方可以用兩個OMEGA相乘得到,但由於容易產生浮點誤差,加上數值固定,因此也直接做成常量
 *  
 */
ComplexConsts.OMEGA_SQUARE = new ComplexNum(1 / 2, Math.sqrt(3) / 2);		

/**
 * OMEGA的3次方恰好等於實數1,就最好不要再自己用三個OMEGA來相乘了
 *  
 */
ComplexConsts.OMEGA_CUBIC = new ComplexNum(1, 0);

/**
 * OMEGA的0次方,純屬為了讓程式碼清晰而定義的常量
 *  
 */
ComplexConsts.OMEGA_ZERO = new ComplexNum(1, 0);

/**
 * 虛數單位i
 *  
 */
ComplexConsts.I = new ComplexNum(0, 1);

AngleUtil.js

function AngleUtil()
{
	
}

AngleUtil.getCosDic = function()
{
	if(AngleUtil._cosDic == undefined)
	{
		AngleUtil._cosDic = [];
		AngleUtil._cosDic[60] = AngleUtil._cosDic[-60] = AngleUtil._cosDic[300] = AngleUtil._cosDic[-300] = 0.5;
		AngleUtil._cosDic[120] = AngleUtil._cosDic[-120] = AngleUtil._cosDic[240] = AngleUtil._cosDic[-240] = -0.5;
		AngleUtil._cosDic[0] = AngleUtil._cosDic[360] = 1;
		AngleUtil._cosDic[180] = AngleUtil._cosDic[-180] = -1;
		AngleUtil._cosDic[90] = AngleUtil._cosDic[270] = AngleUtil._cosDic[-90] = AngleUtil._cosDic[-270] = 0;	
	}
	return AngleUtil._cosDic;
}

AngleUtil.getSinDic = function()
{
	if(AngleUtil._sinDic == undefined)
	{
		AngleUtil._sinDic = [];
		AngleUtil._sinDic[30] = AngleUtil._sinDic[150] = AngleUtil._sinDic[-330] = AngleUtil._sinDic[-210] = 0.5;
		AngleUtil._sinDic[-30] = AngleUtil._sinDic[-150] = AngleUtil._sinDic[210] = AngleUtil._sinDic[330] = -0.5;
		AngleUtil._sinDic[90] = AngleUtil._sinDic[-270] = 1;
		
		AngleUtil._sinDic[270] = AngleUtil._sinDic[-90] = -1;
		AngleUtil._sinDic[180] = AngleUtil._sinDic[-180] = AngleUtil._sinDic[0] = AngleUtil._sinDic[360] = 0;
	}
	return AngleUtil._sinDic;
}

真沒想到,弄完這些輔助類就把文章撐那麼長了。那我就分開兩篇寫好了。

其實我可以弄附件,但有些東西沒測好,附件又似乎不太好改,那就還是直接貼吧。

下篇我貼出解方程的核心類。