[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation

Introduction

** Interpolate using scipy's interpolate lagrange method. ** **

** Interpolate a given N + 1 dataset ($ x_i $, $ y_i $) (i = 0,2, 3, ..., N) by Nth order polynomial ([Lagrange interpolation](https:: //ja.wikipedia.org/wiki/%E3%83%A9%E3%82%B0%E3%83%A9%E3%83%B3%E3%82%B8%E3%83%A5%E8%A3 % 9C% E9% 96% 93)) Implement the program in Python3. ** **

Consider $ y = 1 / (1 + x ^ 2) $ as an example. Sample an 11-point dataset ($ x_i $, $ y_i $) and interpolate it with a 10th-order polynomial. This function does not work with Lagrange polynomial [Runge phenomenon](https://ja.wikipedia.org/wiki/%E3%83%AB%E3%83%B3%E3%82%B2%E7%8F%BE%E8 It is known as an example of causing% B1% A1).


Contents

Code 1: Use scipy / numpy. Minimal code (** Click here if you are in a hurry **) Code 2: Calculate the coefficients for interpolation in the code. Click here to study the method.


Calculation code (1): Use scipy / numpy

code(1)scipyを利用してラクをするcode。


from scipy.interpolate import lagrange
import numpy as np
import matplotlib.pyplot as plt

##Main
x =np.linspace(-5,5,num=11) #[-5,5]Divide the range of into 11 equal parts and store in x
y = 1.0/(1.0+x**2)          #The function considered in this example, y= 1/(1+x^2)Define
f_Lag=lagrange(x,y) #scipy.interpolate.Lagrange interpolation execution by lagrange
##

#for plot
xnew =np.linspace(-5,5,num=51) # [-5,5]Divide the range of into 51 equal parts and store in xnew
plt.plot(x, y, 'o',  xnew, f_Lag(xnew), '-')   #Raw data"o"So, Lagrange interpolated line('-')Draw with.
plt.legend(['Raw data','Lagrange'], loc='best')  #specification of legend
plt.xlim([-6, 6])  #x-axis plot range
plt.ylim([0, 1.4]) #y-axis plot range
plt.show()

Result (1)

t.png

11 data points sampled by blue. The orange line is Lagrange polynomial.


Calculation code (2) Calculate the coefficient for interpolation with python

code(2)code中に直接数値計算を実行する



"""
interpolation: ラグランジュinterpolation
example: y = 1/(1+x**2):Section-11 points from 5 to 5 are sampled and Lagrange interpolated.
"""
from math import pi,e, log, factorial
import matplotlib.pyplot as plt


### 
def g(i, x): #Calculation of coefficients of Lagrange interpolation. Main.
    dum=1.0
    for j in range(len(x_lis)):
        if j != i :
            dum *= (x-x_lis[j])/(x_lis[i]-x_lis[j])
    return dum

#
def fLag(x,m): #Lagrange interpolation
    dum=0.0
    for j in range(m):
        dum += y_lis[j]*g(j, x)
    
    return dum
###
    
##Dataset construction for the example
m = 11  #  x= -11-point sampling at regular intervals from 5 to 5
x_lis = []
y_lis = []
def yy(x):
    return 1/(1+x**2) #Example function y= 1/(1+x**2)

for k in range(m):
    xm = -5.0 + 10.0*k / m
    x_lis.append(xm)
    y_lis.append(yy(xm))
    
plt.plot(x_lis,y_lis, 'o',label='Row data')
##

###Lagrange interpolation execution
mm = 5000
y_Laglis = []
xx_lis = []
for k  in range(mm):
    xm = -5.0 + 10.0*k / mm
    xx_lis.append(xm)

y_lis_exact=[]
for j in range(mm):
    y_Laglis.append(fLag(xx_lis[j],m))
    y_lis_exact.append(yy(xx_lis[j]))

    
#plot
plt.grid(True)
plt.xlabel('x',fontsize=24)
plt.ylabel('f(x)',fontsize=24)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)
plt.plot(xx_lis,y_Laglis, color='Red',label='Lagrange')
plt.plot(xx_lis,y_lis_exact, color='Black',label='Exact')
plt.legend(loc='upper left')
plt.show()¥

Result (2)

Lag.png

11 data points sampled by blue. The red line is Lagrange polynomial. The black line represents the exact solution.

** As the order (data points) is increased, the approximation of the central part becomes better, but the approximation of the part near both ends becomes worse. It is known as an example in which the error increases as the order is increased and finally diverges to infinity. This is because equinox equinoxes are used, and it is known that it can be avoided if the coordination of the equinoxes is rough near the center and finely divided at both ends [1]. ** **


Addendum: Conditions for evenly spaced Lagrange interpolation to approach the original function f (x) as the score increases.

** For interval [-1,1]: When the domain of $ f (x) $ is analyzed and connected to the complex plane, its singularity $ z $ is **

|(1+z)^{1+z}| \leqq 4 e^{\pi |Im(z)|}

** Does not meet **. In this example, the singularity z of $ f (z) = 1 / (1 + z ^ 2) $ is $ ± i $. The right side of the inequality is $ 4e ^ \ pi = 92.562770531 ... , The left side is|(1+i)^{1+i}| = 0.644793...$ Therefore, since the above inequality is ** satisfied **, the error of Lagrange interpolation (polynomial interpolation) using evenly spaced equinoxes increases as the number of equinoxes increases.

If f (x) is a continuous function in the interval [a, b], the interpolation polynomial approaches f (x) at the limit where N is infinite if the equinox $ x_n $ is selected as follows [2]. ..

x_n = (a+b)/2 + (b-a)(sin(-\pi/2+n\pi/N))/2, $(n = 0, 1, 2, ...,N) $

In addition, polynomial interpolation that selects the zero point of an orthogonal polynomial (for example, Chebisif polynomial) as the equinox does not cause the Rung phenomenon as seen above [1].

References

[1] Masatake Mori, ["Numerical Analysis Second Edition"](https://www.amazon.co.jp/%E6%95%B0%E5%80%A4%E8%A7%A3%E6%9E % 90-% E5% 85% B1% E7% AB% 8B% E6% 95% B0% E5% AD% A6% E8% AC% 9B% E5% BA% A7-% E6% A3% AE-% E6% AD% A3% E6% AD% A6 / dp / 4320017013), Kyoritsu Shuppan, 2002. [2] Masao Iri and Kazutake Fujino, ["Common knowledge of numerical calculation"](https://www.amazon.co.jp/%E6%95%B0%E5%80%A4%E8%A8%88 % E7% AE% 97% E3% 81% AE% E5% B8% B8% E8% AD% 98-% E4% BC% 8A% E7% 90% 86-% E6% AD% A3% E5% A4% AB / dp / 4320013433), Kyoritsu Shuppan, 1985.

Recommended Posts

[Scientific / technical calculation by Python] Lagrange interpolation, numerical calculation
[Scientific / technical calculation by Python] Sum calculation, numerical calculation
[Scientific / technical calculation by Python] 3rd order spline interpolation, scipy
[Scientific / technical calculation by Python] Monte Carlo integration, numerical calculation, numpy
[Scientific / technical calculation by Python] Numerical integration, trapezoidal / Simpson law, numerical calculation, scipy
[Scientific / technical calculation by Python] Solving simultaneous linear equations, numerical calculation, numpy
[Scientific / technical calculation by Python] 2D random walk (drunk walk problem), numerical calculation
[Scientific / technical calculation by Python] Inverse matrix calculation, numpy
[Scientific / technical calculation by Python] Histogram, visualization, matplotlib
[Scientific / technical calculation by Python] Logarithmic graph, visualization, matplotlib
[Scientific / technical calculation by Python] Polar graph, visualization, matplotlib
[Scientific / technical calculation by Python] Solving second-order ordinary differential equations by Numerov method, numerical calculation
[Scientific / technical calculation by Python] Numerical calculation to find the value of derivative (differential)
[Scientific / technical calculation by Python] Basic operation of arrays, numpy
[Scientific / technical calculation by Python] Vector field visualization example, electrostatic field, matplotlib
[Scientific / technical calculation by Python] Numerical solution of second-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] List of matrices that appear in Hinpan in numerical linear algebra
[Scientific / technical calculation by Python] Numerical solution of one-dimensional harmonic oscillator problem by velocity Verlet method
[Scientific / technical calculation by Python] Numerical solution of eigenvalue problem of matrix by power method, numerical linear algebra
[Scientific / technical calculation by Python] Fitting by nonlinear function, equation of state, scipy
[Scientific / technical calculation by Python] Calculation of matrix product by @ operator, python3.5 or later, numpy
[Scientific / technical calculation by Python] Solving ordinary differential equations, mathematical formulas, sympy
Lagrange interpolation in python
Numerical calculation with Python
[Scientific / technical calculation by Python] Solving (generalized) eigenvalue problem using numpy / scipy, using library
[Scientific / technical calculation by Python] Drawing animation of parabolic motion with locus, matplotlib
[Scientific / technical calculation by Python] Analytical solution to find the solution of equation sympy
[Scientific / technical calculation by Python] Plot, visualize, matplotlib 2D data with error bars
[Scientific / technical calculation by Python] Drawing of 3D curved surface, surface, wireframe, visualization, matplotlib
[Scientific / technical calculation by Python] Solving one-dimensional Newton equation by the 4th-order Runge-Kutta method
[Scientific / technical calculation by Python] Solving the boundary value problem of ordinary differential equations in matrix format, numerical calculation
[Scientific / technical calculation by Python] Plot, visualization, matplotlib of 2D data read from file
[Scientific / technical calculation by Python] Drawing, visualization, matplotlib of 2D (color) contour lines, etc.
[Scientific / technical calculation by Python] Numerical solution of one-dimensional and two-dimensional wave equations by FTCS method (explicit method), hyperbolic partial differential equations
[Science / technical calculation by Python] Numerical solution of first-order ordinary differential equations, initial value problem, numerical calculation
[Scientific / technical calculation by Python] List of usage of (special) functions used in physics by using scipy
[Scientific and technical calculation by Python] Drawing of fractal figures [Sierpinski triangle, Bernsley fern, fractal tree]
[Scientific / technical calculation by Python] Wave "beat" and group velocity, wave superposition, visualization, high school physics
[Scientific / technical calculation by Python] Monte Carlo simulation of thermodynamics of 2D Ising spin system by Metropolis method
Scientific / technical calculation by Python] Drawing and visualization of 3D isosurface and its cross section using mayavi
Create a new Python numerical calculation project
[Scientific / technical calculation by Python] Derivation of analytical solutions for quadratic and cubic equations, mathematical formulas, sympy
[Scientific / technical calculation by Python] Numerical solution of two-dimensional Laplace-Poisson equation by Jacobi method for electrostatic potential, elliptic partial differential equation, boundary value problem
[Scientific / technical calculation by Python] Numerical solution of one-dimensional unsteady heat conduction equation by Crank-Nicholson method (implicit method) and FTCS method (positive solution method), parabolic partial differential equation
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (1), well-type potential, quantum mechanics
Video frame interpolation by deep learning Part1 [Python]
Calculation of technical indicators by TA-Lib and pandas
[Scientific / technical calculation by Python] Generation of non-uniform random numbers giving a given probability density function, Monte Carlo simulation
[Scientific / technical calculation by Python] Solving one-dimensional Schrodinger equation in stationary state by shooting method (2), harmonic oscillator potential, quantum mechanics
[Python] Numerical calculation of two-dimensional thermal diffusion equation by ADI method (alternate direction implicit method)
Numerical calculation of compressible fluid by finite volume method
Start numerical calculation in Python (with Homebrew and pip)