Нужна помощь в построении этого кода кусочно-линейной интерполяции в Python - PullRequest
Новая
       89

Нужна помощь в построении этого кода кусочно-линейной интерполяции в Python

0 голосов
/ 19 июня 2020

Я пробую интерполировать данные с помощью простой линейной функции, интерполяционного полинома Лагранжа. Мне удалось получить необходимые уравнения, однако я не могу понять, как их построить по частям. Я понимаю, что использование sympy - не лучший путь вперед, но я новичок и хотел увидеть, как выглядят мои уравнения.

Как я могу построить график в matplotlib без необходимости вручную вводить уравнения в конце? 

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

x = sym.Symbol('x')

year = np.arange(1960,2020,10)
pop = [179323,203302,226542,249633,281422,308746]

def lgn(a,b): #this ideally should be taking a value of x where you'd like to interpolate.
    result = []
    for i in range(1,len(a)):
        L0 = (x - a[i])/(a[i-1] - a[i])
        L1 = (x - a[i-1])/(a[i] - a[i-1])
        temp = (L0 * b[i-1]) + (L1 * b[i])
        result.append(temp)
    return result

lgn(year,pop) #result gives a list of linear equations between each year value. 

[23979*x/10 - 4520561,
 2324*x - 4374978,
 23091*x/10 - 4345476,
 31789*x/10 - 6076378,
 13662*x/5 - 5183378]

#plotting for each interval. this is what I am trying to code.  
x1 = np.linspace(year[0],year[1],10)
y1 = 23979 * x1/10 - 4520561

x2 = np.linspace(year[1],year[2],10)
y2 = 2324*x2 - 4374978

x3 = np.linspace(year[2],year[3],10)
y3 = 23091*x3/10 - 4345476

x4 = np.linspace(year[3],year[4],10)
y4 = 31789*x4/10 - 6076378

x5 = np.linspace(year[4],year[5],10)
y5 = 13662*x5/5 - 5183378

plt.plot(year,pop,'ro',x1,y1,x2,y2,x3,y3,x4,y4,x5,y5 )

1 Ответ

1 голос
/ 19 июня 2020

Для преобразования симпозиумного выражения в числовое значение с помощью expr.subs(x, 123).evalf(). Обратите внимание, что это работает только для замены одного значения. Для работы с массивами sym.lambdify() может преобразовать выражение в функцию, которая понимает numpy массивы, которые затем можно использовать для построения графиков.

Вот пример кода: 

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

def lgn(a, b):
    result = []
    for i in range(1, len(a)):
        L0 = (x - a[i]) / (a[i - 1] - a[i])
        L1 = (x - a[i - 1]) / (a[i] - a[i - 1])
        temp = (L0 * b[i - 1]) + (L1 * b[i])
        result.append(temp)
    return result

x = sym.Symbol('x')

year = np.arange(1960, 2020, 10)
pop = [179323, 203302, 226542, 249633, 281422, 308746]

equations = lgn(year, pop)
for i in range(1, len(year)):
    xi = np.linspace(year[i - 1], year[i], 10)
    yi_np = sym.lambdify(x, equations[i - 1])
    yi = yi_np(xi)
    plt.plot(xi, yi)
plt.plot(year, pop, 'ro')
plt.show()

example plot

Вот подход с использованием sympy's Piecewise: 

import numpy as np
import matplotlib.pyplot as plt
import sympy as sym

def line_2points(a0, b0, a1, b1):
    L0 = (x - a1) / (a0 - a1)
    L1 = (x - a0) / (a1 - a0)
    return L0 * b0 + L1 * b1

x = sym.Symbol('x')

year = np.arange(1960, 2020, 10)
pop = [179323, 203302, 226542, 249633, 281422, 308746]

eq = 0
for i in range(1, len(year)):
    eq = sym.Piecewise( (line_2points(year[i-1], pop[i-1], year[i], pop[i]), (x >= year[i-1] ) &  (x <= year[i] ) ),
                        (eq, True) )
# sym.plot(eq, (x, year[0], year[-1])) # this also works, but the visualization is much harder to customize

eq_np = sym.lambdify(x, eq)
xs = np.linspace(year[0], year[-1], 200)
plt.plot(xs, eq_np(xs))
plt.plot(year, pop, 'ro')
plt.show()
...