scipy 使用lambda函数定义显示表达式

sf6xfgos  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(197)

我使用下面的代码来求分段多项式函数

import scipy
from scipy.interpolate import UnivariateSpline, splrep

x = np.array([0., 0.75, 1.8, 2.25, 3.75, 4.5, 6.45, 6.75, 
          7.5, 8.325, 10.875, 11.25, 12.525, 12.75, 15., 
          20.85, 21.])
y = np.array([2.83811035, 2.81541896, 3.14311655, 
          3.22373554, 3.43033456, 3.50433385, 
          3.66794514, 3.462296,   3.59480959, 
          3.56250726, 3.6209845,  3.63034523, 
          3.68238915, 3.69096892, 3.75560395, 
          3.83545191, 3.90419498])

k = 3  # polynomial order
spl = UnivariateSpline(x, y, k=3, s=0.09)
xs = np.linspace(x.min(), x.max(), 100)
plt.plot(x, y, 'ro', ms=5)
plt.plot(xs, spl(xs), 'cyan', lw=5, alpha=0.3)

# get spline coeffs and knots

tck = (spl._data[8], spl._data[9], k)  # tck = (knots, coefficients, degree)
p = scipy.interpolate.PPoly.from_spline(tck)

# plot each segment and return knots and coeffs

for idx, i in enumerate(range(k, len(spl.get_knots()) + k - 1)):
    xs = np.linspace(p.x[i], p.x[i + 1], 100)
    plt.plot(xs, np.polyval(p.c[:, i], xs - p.x[i]))
    print("knot ", p.x[i], " to ", p.x[i + 1])
    print("coeffs ", p.c[:, i], "\n")
    f0 = lambda x: p.c[0, i] * (x - p.x[i])**3 + p.c[1, i] * (x - p.x[i])**2 + p.c[2, i] * (x - p.x[i]) + p.c[3, i]
    f0 = lambda x: [p.c[:, i] * (x - p.x[i])**(3 - i) for i in range(k + 1)]  # k = degree i.e no. of. coeffs = degree +1
    print(f0)
plt.show()

print(f0)输出
相反,我希望显示多项式表达式。
关于如何做到这一点的建议将会非常有帮助。
编辑:
我在一个帖子https://stackoverflow.com/a/60114991/8281509中找到了以下内容。但不知道如何将其用于我的情况。

from sympy import lambdify, bspline_basis_set
from sympy.abc import u

basis = bspline_basis_set(tck[2], tck[0],  u)
for i, b in enumerate(basis):
    print(f"Basis {i} :", b)

返回16个基函数,而我只期望3个

62lalag4

62lalag41#

您可以使用如下的字串格式:

for idx, i in enumerate(range(k, len(spl.get_knots()) + k - 1)):
    xs = np.linspace(p.x[i], p.x[i + 1], 100)
    plt.plot(xs, np.polyval(p.c[:, i], xs - p.x[i]))
    print("knot ", p.x[i], " to ", p.x[i + 1])
    print("coeffs ", p.c[:, i], "\n")
    print(f"{p.c[0, i]} * {(x - p.x[i])}**3 + {p.c[1, i]} * {(x - p.x[i])}**2 + {p.c[2, i]} * {(x - p.x[i])} + {p.c[3, i]}")

输出量:

0.00041014998689577995 * [ 0.     0.75   1.8    2.25   3.75   4.5    6.45   6.75   7.5    8.325
 10.875 11.25  12.525 12.75  15.    20.85  21.   ]**3 + -0.015832105893445334 * [ 0.     0.75   1.8    2.25   3.75   4.5    6.45   6.75   7.5    8.325
 10.875 11.25  12.525 12.75  15.    20.85  21.   ]**2 + 0.20309856945144328 * [ 0.     0.75   1.8    2.25   3.75   4.5    6.45   6.75   7.5    8.325
 10.875 11.25  12.525 12.75  15.    20.85  21.   ] + 2.8044698839070015

数组的内容会被打印出来。如果你不需要它,就用文字u来替换它,如下所示:

print(f"{p.c[0, i]} * u**3 + {p.c[1, i]} * u**2 + {p.c[2, i]} * u + {p.c[3, i]}")

输出量:

0.00041014998689577995 * u**3 + -0.015832105893445334 * u**2 + 0.20309856945144328 * u + 2.8044698839070015

相关问题