︠b96482d5-8c81-4c9e-81cc-b56888bfcd6fi︠
%hide
# by Mao Ziyang (mao-ziyang@yeah.net)
var('x')
@interact
def interpolate(pts = input_box(default=[[-2,2],[0,1],[1,2],[2,7]], label='Points'),
interpolate_rule=selector(['Lagrange', 'Newton'], nrows=1, label='Interpolation rule'),
show_basis=checkbox(False, "Show basis functions")):
n = len(pts)-1 #插值多项式次数
if interpolate_rule == 'Lagrange':
l = [] #插值基函数
for i in range(n+1):
xvalues = [pt[0] for pt in pts]
deleted_x = xvalues.pop(i)
numerator = prod((x - j) for j in xvalues)
denominator = numerator(x = deleted_x)
l.append(numerator / denominator)
p(x) = sum([l[i] * pts[i][1] for i in range(n+1)])
html(r'''
Lagrange Interpolation
\begin{align*}
%s = & %s\\
= & %s
\end{align*}
''' % ('L_%s(x)'%n,latex(p(x)),latex(expand(p(x))))
)
elif interpolate_rule =='Newton':
l = [] #插值基函数
a = [pt[1] for pt in pts]
xvalues = [pt[0] for pt in pts]
for k in range(1, n+1):
for i in range(n, k-1, -1):
a[i] = (a[i] - a[i-1])/(xvalues[i] - xvalues[i-k])
p(x) = a[0]
l.append(1)
for i in range(1, n+1):
l.append(l[i-1]*(x-xvalues[i-1]))
p(x) = p(x) + a[i]*l[i]
#另一种计算方法，计算量更小
#for k in range(1,n+1):
# for i in range(k, n+1):
# a[i] = (a[i] - a[k-1])/(xvalues[i] - xvalues[k-1])
#p(x) = a[n]
#l.append(1)
#for k in range(1,n+1):
# p(x) = a[n-k] + (x - xvalues[n-k]) * p(x)
# l.append((x-xvalues[n-k])*l[k-1])
html(r'''
Newton Interpolation
\begin{align*}
%s = & %s\\
= & %s
\end{align*}
''' % ('N_%s(x)'%n,latex(p(x)),latex(expand(p(x))))
)
xvalues = [pt[0] for pt in pts]
fig = scatter_plot(pts)+plot(p(x), [min(xvalues),max(xvalues)])
fig.show()
if show_basis:
fig2 = plot([l[i] for i in range(n+1)], [min(xvalues), max(xvalues)], fill=true)
if interpolate_rule == 'Lagrange':
fig2 += scatter_plot([(pt[0],1) for pt in pts])
fig2.show()
︡b5873ba1-e22a-47e2-b8d7-7f2b15fbacc6︡{"html": "\n "}︡