【问题标题】:Taylor Expansion for given dataset without functional form没有函数形式的给定数据集的泰勒展开
【发布时间】:2015-10-02 23:38:32
【问题描述】:

我有 (x,y) 数据集,它是连续且可微的。确切的函​​数形式是未知的。我想在某个时候扩展图表。我试过使用algopy/Adipy。问题是他们需要功能形式。

我附上algopy的示例代码。

 import numpy; from numpy import sin,cos
 from algopy import UTPM

 def f(x):
     return sin(cos(x) + sin(x))

D = 100; P = 1
x = UTPM(numpy.zeros((D,P)))
x.data[0,0] = 0.3
x.data[1,0] = 1
y = f(x)
print('coefficients of y =', y.data[:,0])

其中 D 是多项式的阶数。

我尝试使用以下(x1 和 y1 是一维数组):

from scipy.interpolate import interp1d
f1 = interp1d(x1, y1, kind='cubic')
def f(x):
    temp1=f1(x)
    return np.float64(temp1)

但是,插值似乎不接受 UTPM 返回的 x 数据类型。

错误信息:

Traceback (most recent call last):
      File "tay.py", line 26, in <module>
        y = f(x)
      File "tay.py", line 15, in f
        temp1=f1(x)
      File "/usr/lib/python2.7/dist-packages/scipy/interpolate/polyint.py", line 54, in __call__
        y = self._evaluate(x)
      File "/usr/lib/python2.7/dist-packages/scipy/interpolate/interpolate.py", line 449, in _evaluate
        y_new = self._call(self, x_new)
      File "/usr/lib/python2.7/dist-packages/scipy/interpolate/interpolate.py", line 441, in _call_spline
        return spleval(self._spline, x_new)
      File "/usr/lib/python2.7/dist-packages/scipy/interpolate/interpolate.py", line 919, in spleval
        res[sl] = _fitpack._bspleval(xx,xj,cvals[sl],k,deriv)
    TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'

【问题讨论】:

标签: python scipy taylor-series


【解决方案1】:

对定义在离散点的数据集进行泰勒展开是没有意义的。特别是下面的介词是错误的,

我有 (x,y) 数据集,它是连续且可微的。确切的函​​数形式未知。

如果你将一些插值过程与你的数据集相关联,你只能有一个连续函数,但这也将修复一般函数形式。

例如,假设我们使用分段三次插值,如问题所示。这意味着泰勒展开式将受到用于插值的三次多项式系数的有效约束(并且最多可以是 3 阶)。此外,另一个插值程序将产生不同的泰勒展开。

一般来说,结果主要取决于插值例程,而不是您的数据。这是因为泰勒展开依赖于函数的局部行为,而该函数不包含在您的 (x,y) 数据集中。

相反,您可以使用某个阶的多项式在本地拟合数据,这将产生与采样数据的泰勒展开式等效的结果。

【讨论】:

    【解决方案2】:

    我一直在寻找相同的东西,所以我实现了这个:

    import numpy as np
    import matplotlib.pyplot as plt
    from math import factorial as f
    
    def dxdy(x,y,order): 
        dy = y
        for k  in range(order+1):
            print(k)
            dx = (x[-1]-x[0])/len(x)
            if k == 1:
                dy = y
            elif k % 2 == 0:
                dy = (dy[1:]-dy[:-1])/dx
                x = x[:-1]       
            elif k % 2 != 0:
                dy = (dy[1:]-dy[:-1])/dx
                x = x[1:]
        return dy
    
    def taylor(x,y,n):
        a = x[int(len(x)/2)+1]
        center = int(len(x)/2)+1o
        #plt.plot(y)
        #plt.ylim(min(y),max(y))
        for k in range(n+1):
            print(k)
            if k == 0:
                y_hat = (y[center]*((x-a)**k))/f(k)
                #plt.plot(y_hat)
            else:
                y_hat += (dxdy(x,y,k+1)[center]*((x-a)**k))/f(k)
                #plt.plot(y_hat)
            #plt.plot(y)
        return y_hat
    
    points = 101
    x = np.linspace(-3*np.pi,3*np.pi,points)
    y = 1/(1+np.exp(-x))
    y = np.cos(x)#*x#(x**4)
    center = int(points/2)
    for k in range(21):
        y_hat = taylor(x,y,k)
        plt.figure(figsize=(8,4))
        plt.ylim(min(y)*1.1,max(y)*1.1)
        plt.xlim(min(x),max(x))
        plt.plot(x,y)
        plt.plot(x,y_hat,c='red')
        plt.legend(['cs(x)','taylor, k= '+str(k)],loc='upper right')
        plt.title('cos(x)') 
        plt.savefig('cos'+str(k)+'.png')
    

    【讨论】:

      猜你喜欢
      • 2014-06-11
      • 2017-06-14
      • 2013-10-09
      • 2017-04-02
      • 2021-07-10
      • 2020-01-05
      • 2020-04-10
      • 1970-01-01
      • 2019-06-15
      相关资源
      最近更新 更多