Find the solution of the nth-order equation in python

What is the nth-order equation?

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.

algorithm

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.

Implementation

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.

References

Proof of relationship between companion matrix and characteristic polynomial

Recommended Posts

Find the solution of the nth-order equation in python
Find the divisor of the value entered in python
Solving the equation of motion in Python (odeint)
Implement the solution of Riccati algebraic equations in Python
[Scientific / technical calculation by Python] Analytical solution to find the solution of equation sympy
Find out the apparent width of a string in python
Have the equation graph of the linear function drawn in Python
Python --Find out number of groups in the regex expression
Find the eigenvalues of a real symmetric matrix in Python
Check the behavior of destructor in Python
The result of installing python in Anaconda
The basics of running NoxPlayer in Python
In search of the fastest FizzBuzz in Python
Double pendulum equation of motion in python
Find the numerical solution of the second-order ordinary differential equation with scipy
Find the maximum Python
[Python] Sort the list of pathlib.Path in natural sort
Get the caller of a function in Python
Match the distribution of each group in Python
View the result of geometry processing in Python
the zen of Python
Make a copy of the list in Python
Find the number of days in a month
List find in Python
The story of reading HSPICE data in Python
[Note] About the role of underscore "_" in Python
About the behavior of Model.get_or_create () of peewee in Python
Output in the form of a python array
Find the geometric mean of n! Using Python
[Python] Find the transposed matrix in a comprehension
python xlwings: Find the cell in the last row
How to find the coefficient of the trendline that passes through the vertices in Python
Experience the good calculation efficiency of vectorization in Python
How to get the number of digits in Python
Find out the location of Python class definition files.
[python] Get the list of classes defined in the module
The story of FileNotFound in Python open () mode ='w'
Learn the design pattern "Chain of Responsibility" in Python
Find the part that is 575 from Wikipedia in Python
Get the size (number of elements) of UnionFind in Python
Not being aware of the contents of the data in python
Find the Hermitian matrix and its eigenvalues in Python
Reproduce the execution example of Chapter 4 of Hajipata in Python
Let's use the open data of "Mamebus" in Python
Implemented the algorithm of "Algorithm Picture Book" in Python3 (Heapsort)
[Python] Outputs all combinations of elements in the list
Get the URL of the HTTP redirect destination in Python
A reminder about the implementation of recommendations in Python
Reproduce the execution example of Chapter 5 of Hajipata in Python
To do the equivalent of Ruby's ObjectSpace._id2ref in Python
Check the asymptotic nature of the probability distribution in Python
Download the file in Python
About the ease of Python
Equivalence of objects in Python
Implementation of quicksort in Python
About the features of Python
Let's find pi in Python
The Power of Pandas: Python
Find the maximum python (improved)
I used Python to find out about the role choices of the 51 "Yachts" in the world.
I tried to find the entropy of the image with python