3

我有 (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'
4

2 回答 2

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')

在此处输入图像描述

于 2020-02-11T13:03:14.720 回答
2

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

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

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

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

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

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

于 2015-07-14T12:19:11.110 回答