0

I am trying my hand on interpolating data using a simple linear function, Lagrange Interpolating Polynomial. I have managed to get the required equations, however, I am not able to figure out how to plot it piece-wise. I do understand using sympy is not the best way forward, but I am a noob and I wanted to see how my equations look.

How can I make it plot in matplotlib without having the need to manually type the equations at the end?

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 Answers1

1

To convert a sympy expression to a numerical value with expr.subs(x, 123).evalf(). Note that this only works for substituting a single value. To work with arrays, sym.lambdify() can convert the expression to a function that understands numpy arrays, which then can be used for plotting.

Here is some example code:

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

Here is an approach using 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()
JohanC
  • 71,591
  • 8
  • 33
  • 66
  • Thanks. I know sympy might now be the best way to do this, what could be other possible solutions in that case? Do you think this plot code will work for a quadratic/cubic spline functions as well? – Harsh Dhiman Jun 19 '20 at 11:04
  • See [Construct a symbolic interpolating spline through given points using SymPy ](https://stackoverflow.com/questions/42602553/construct-a-symbolic-interpolating-spline-through-given-points-using-sympy) about a way to construct an interpolating spline with sympy. You can lambdify such a spline similarly as here. – JohanC Jun 19 '20 at 11:19
  • Thanks, that was very helpful. – Harsh Dhiman Jun 19 '20 at 13:37