x^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 =0
I came across a situation where I wanted to solve. I feel that this is generally called the nth-order equation. (Is it also called a one-variable algebraic equation?) It is known that there is no algebraic solution for those of order 5 and above. However, it seems that the fifth order can be expressed by using various special functions, and the approximate solution can be calculated numerically. It seems that it is normal to use numerical calculation for the 4th order and above. It seems that Matlab has a function called roots that performs numerical calculation of algebraic equations, but I could not find a library in python. That's why I decided to implement it.
When it comes to finding a solution, Newton's method seems to come up first, but it seems that the DKA method is known as an algorithm similar to Newton's method for finding solutions to nth-order equations. However, implementing this algorithm is tedious, and if you look at Implementation of Matlab, you can use a method that results in eigenvalue calculation. There seems to be. I decided to use this this time. When calculating eigenvalues according to the definition
\mathrm{det}(sI-A)=0
However, this is an nth-order equation when the matrix is nth-order. There seem to be various methods for calculating eigenvalues, and numpy also has a function called np.linalg.eig. Therefore, it is only necessary to know the matrix that has the equation for which the solution is to be obtained as the characteristic polynomial. For this, the following companion matrix is known.
A=
\left(
\begin{array}{ccccc}
0 & 1 & 0&\ldots&0 \\
\vdots & \ddots & 1&\ddots&\vdots \\
\vdots & \ldots & \ddots&\ddots&0\\
0&\ldots&\ldots&0&1\\
-a_0&-a_1&\ldots&-a_{n-2}&-a_{n-1}
\end{array}
\right)
The eigen equation of this is
|sI-A| = s^n+a_{n-1}s^{n-1}+a_{n-2}s^{n-2}+\ldots+a_1s+a_0
Will be. Therefore, if all the eigenvalues of the above companion matrix are calculated for the nth-order equation with arbitrary coefficients, the solution can be obtained. Please refer to here for the DKA method introduced first. Implementation was also down.
n-poly.py
import numpy as np
def solve(vec,is_complex=False):
dim =len(vec)
if is_complex:
A = np.zeros((dim,dim),dtype=complex)
else:
A = np.zeros((dim,dim))
A[np.arange(dim-1),1+np.arange(dim-1)] =1
A[-1,:] = -vec
ans,vec = np.linalg.eig(A)
return ans
vec0 =np.array([-120,274,-225,85,-15])
vec1 =np.array([1,0])
vec2 =np.array([-1,5,-10,10,-5])
print(solve(vec0))
print(solve(vec1))
print(solve(vec2))
[ 5. 4. 3. 2. 1.]
[ 0.+1.j 0.-1.j]
[ 1.00079742+0.00057977j 1.00079742-0.00057977j 0.99969502+0.00093688j
0.99969502-0.00093688j 0.99901512+0.j ]
Note that if the coefficient contains an imaginary number, the coefficient will be converted to a real number without setting is_complex to True. If dtype = complex is always set, it seems that the real eigenvalue is difficult to calculate as a real eigenvalue, so it is usually calculated by float (Is the algorithm different depending on the dtype?)
Each example
(x-1)(x-2)(x-3)(x-4)(x-5)=x^5-15x^4+85x^3-225x^2+274x-120=0\\
x^2+1=(x-i)(x+i)=0\\
(x-1)^5=0
Is being calculated. It can be seen that even quintic equations without a solution formula are calculated with reasonable accuracy, including multiple solutions.
Proof of relationship between companion matrix and characteristic polynomial
Recommended Posts