The integral of the complex trapezoid--numerical calculation with Python

There are many ways to use the program to find the integral, this article is mainly about Newton-Cortez formula.

The easiest thing for a classmate to learn about interpolation is to use an interpolation function instead of an integral function to find an integral, but in reality it doesn't work in most scenarios.

The interpolation function is usually a polynomial of not more than n times, if the interpolation function to find the integral, it will introduce the problem of the higher order polynomial to seek integral. This will take the original quadrature problem to another quadrature problem: How to find the integral of the n-th polynomial, and when the number of times become high, there will be a dragon elegy phenomenon, the error may increase, and the high-order interpolation formula may become unstable: detailed reasons do not repeat.

Newton-Cortez Formula the solution to this problem is to divide the large interpolation interval into a small number of interpolation intervals, so that the frequency of the polynomial is not too high. Then by introducing the parametric function

The value range of an item with a power is fixed in a fixed range so that the quadrature of the part with the power of the polynomial becomes a fixed constant, which can only be calculated by hand. This constant can be directly brought into the polynomial quadrature function.

The quadrature interval for x in the above equation is [A, b],h = (b-a)/n, so that the integration interval becomes [0, N], it should be noted from this formula that a large interval is divided into a small interval of n equal length. This section is detailed in any book on numerical calculations!

N is a pre-determined value.

Also, because a large interpolation interval needs to be divided into several small intervals of equal length, and the interpolation and integration are performed on each of these cells, the Newton-Cortez formula at this time is called: the complex Newton-Cortez formula.

And for the different values of N Newton-Cortez have different names: when N=1, called the re-trapezoid formula, the re-trapezoid formula is to see each cell as a trapezoid (h, the bottom is f (t), the bottom is F (t+1)). This is the same as the essence of the integral: infinite separation .

When n=2, the complex Newton-Cortez formula is known as the complex Simpson formula (the one that is not known to the American legal profession, The Simpsons).

In this article, I have implemented a compound trapezoid formula:

　　　　

First write a function to calculate the value of the node function sum that part:

　　　　

"" "@brief: Sum ∑f (XK): XK represents the K-node of the equidistant node, excluding the endpoint        XK = a + kh (k = 0, 1, 2, ...)        The integral interval is [A, b]        @param:      The x-coordinate set (excluding endpoints) of the xk integral interval @param:    func quadrature function @return: The return value is set and "" "Def sum_fun_xk (XK, func) :    return sum ([func (each) for each in XK])

　　

Then you can write the entire function of the quadrature:

"" "@brief: Seek the Func integral:        @param: A  integral interval left endpoint @param:b  integral interval Right endpoint @param:n integral is  divided into n equal parts (complex trapezoidal quadrature requirements) @param: Func  quadrature function @return: Integral value "" "   def Integral (A, B, N, func):    h = (b-a)/float (n)    xk = [A + i*h for I in range (1, N)]    Return H/2 * (func (a) + 2 * SUM_FUN_XK (XK, func) + func (b))

Quite simple.

　　

Test:

When dividing a large interval into two cells:

　　　　

When divided into 20 community rooms:

　　

The integral value of these colored trapezoidal areas.

Test code:

if __name__ = = "__main__":        func = lambda x:x**2    A, B = 2, 8    n =    print integral (A, B, N, func)        " Paint "    Import matplotlib.pyplot as Plt    plt.figure (" Play ")    Ax1 = Plt.subplot (111)    Plt.sca (AX1)        TMPX = [2 + float (8-2)/50 * each for each in range (50+1)]     plt.plot (tmpx, [func (All) for each in tmpx], LineStyle = '-', color= ' black ') for        rang in range (n):        tmpx = [A + float (8-2)/n * rang, A + float (8-2)/n * rang, A + float (8- 2)/n * (rang+1), A + float (8-2)/n * (rang+1)]        tmpy = [0, func (tmpx[1]), func (tmpx[2]), 0]         c = [' R ', ' Y ', ' B ', ' G ']        Plt.fill (tmpx, Tmpy, color=c[rang%4])    Plt.grid (True)    plt.show ()

Note that the n in the above code is not the N in the formula mentioned above, and the first mentioned n refers to inserting n nodes equidistant from each specific interpolation interval (that is, the cell), and the N of the complex trapezoid formula is fixed at 1.

The N in the code divides a large interval into n small intervals.

The integral of the complex trapezoid--numerical calculation with Python