Python数值导数与积分奇异结果

Question

我写了一个脚本来测试一些微分方法。我用正弦波作为测试数据。你知道吗

以下是脚本：

import sys, getopt
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate
from scipy import stats

def d1(x,y):
  dydx=[0]
  dydx.extend(np.diff(y)/np.diff(x))
  xt=zip(x,x[1:],x[2:])
  yt=zip(dydx,dydx[1:],dydx[2:])
  n=0
  dydx=[]
  for i,j in zip(xt,yt):
    X1=(i[0]+i[1])/2
    X2=(i[1]+i[2])/2
    m=(j[2]-j[1])/(X2-X1)
    C=j[1]-m*(i[1]+i[0])/2
    if n==0:
      dydx.append(m*i[0]+C)
    dydx.append(m*i[1]+C)
    if n==len(xt)-1:
      dydx.append(m*i[2]+C)
    n+=1
  return dydx

def d2(x,y):
  x0=x[:-2]
  x1=x[1:-1]
  x2=x[2:]
  y0=y[:-2]
  y1=y[1:-1]
  y2=y[2:]
  f=(x2-x1)/(x2-x0)
  return (1-f)*(y2-y1)/(x2-x1)+f*(y1-y0)/(x1-x0)

def d3(x,y):

  z1 = np.hstack((y[0],y[:-1]))
  z2 = np.hstack((y[1:],y[-1]))
  dx1 = np.hstack((0,np.diff(x)))
  dx2 = np.hstack((np.diff(x),0))
  d=(z2-z1)/(dx2+dx1)
  return d

def d4(x,y):
  dx=np.gradient(x)
  dy=np.gradient(y)
  return dy/dx


def main(argv):
  ifile="/dev/stdin"
  ofile="/dev/stdout"
  plot=False
  graph=False
  snr=False
  ydata=False
  n=1
  dmethod=2

  try:
      opts, args=getopt.getopt(argv,"hi:o:ps:qf:yn:d:r",["ifile=","ofile="])
  except getopt.GetoptError:
    print(sys.argv[0] + ' [-i <inputfile>] [-o <outputfile>] [-d <int> ] [-n <int>] [-p] [-s <imagefile>] [-q] [-f <int>] [-y] [-r]')
    sys,exit(1)
  for opt, arg in opts:
    if opt=='-h':
      print(sys.argv[0] + ' [-i <inputfile>] [-o <outputfile>] [-d <int> ] [-n <int>] [-p] [-s <imagefile>] [-q] [-y]')
      print('\t-n  order of derivative (default=1)')
      print('\t-d  1. modified central difference, 2. central difference (default=2)')
      print('\t-p  display plot')
      print('\t-s  save plot as imagefile (eps, jpeg, jpg, pdf, pgf, png, ps, raw, rgba, svg, svgz, tif, tiff)')
      print('\t-i  inputfile (default=/dev/stdin)')
      print('\t-o  outputfile (default=/dev/stdout)')
      print('\t-q  quiet; suppress output')
      print('\t-r  include signal-to-noise ratio in plot')
      print('\t-y  include y-data in output')
      sys.exit()
    elif opt in ("-i", "--ifile"):
      ifile=arg
    elif opt in ("-o", "--ofile"):
      ofile=arg
    elif opt in ("-q"):
      ofile="/dev/null"
    elif opt in ("-d"):
      dmethod=int(arg)
    elif opt in ("-w"):
      w=int(arg)
    elif opt in ("-p"):
      plot=True
    elif opt in ("-r"):
      snr=True
    elif opt in ("-s"):
      graph=True
      imagefile=arg
    elif opt in ("-y"):
      ydata=True


  data=np.loadtxt(ifile)
  x=data[:,0]
  y=data[:,1]
  yprime=y

  if dmethod==1:
    yprime=d1(data[:,0],data[:,1])
  elif dmethod==2:
    yprime=d3(data[:,0],data[:,1])

# Signal-to-noise calculations
  stonrlist=[]
  stonrlist.append(stats.signaltonoise(y))
  stonrlist.append(stats.signaltonoise(yprime))

# Collect the data and write it
  if ydata:
      outdata=np.array([data[:,0],data[:,1],yprime])
      fmt=['%4.1f','%12.5e','%12.5e']
  else:
      outdata=np.array([data[:,0],yprime])
      fmt=['%4.1f','%12.5e']
  outdata=outdata.T
  np.savetxt(ofile,outdata,fmt)

# Build the plot

  fig=plt.figure()
  fig.set_size_inches(16,9)
  if snr:
    plt.plot(data[:,0],data[:,1],'-',color='black',label=ifile+'\n'+str(stonrlist[0]))
    plt.plot(data[:,0],yprime,'-',color='red',label=ofile+'\n'+str(stonrlist[1]))
  else:
    plt.plot(data[:,0],data[:,1],'-',color='black',label=ifile)
    plt.plot(data[:,0],yprime,'-',color='red',label=ofile)
  plt.legend()
  plt.minorticks_on()
  plt.grid(which='major',linestyle='-',linewidth='0.3')
  plt.grid(which='minor',linestyle=':',linewidth='0.3')
  if plot:
    plt.show()
  if graph:
    fig.savefig(imagefile,dpi=120)

main(sys.argv[1:])

结果如下：

注意振幅有多小。你知道吗

为了测试微分器，我编写了一个积分器。如果他们是任何好的，我应该回到原来的功能（在这种情况下，正弦波）。代码如下：

import sys, getopt
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate

def main(argv):
  ifile="/dev/stdin"
  ofile="/dev/stdout"
  plot=False
  graph=False
  ydata=False

  try:
      opts, args=getopt.getopt(argv,"hi:o:ps:qf:y",["ifile=","ofile="])
  except getopt.GetoptError:
    print(sys.argv[0] + ' [-i <inputfile>] [-o <outputfile>] [-p] [-s <imagefile>] [-q] [-f <int>] [-y]')
    sys,exit(1)
  for opt, arg in opts:
    if opt=='-h':
      print(sys.argv[0] + ' [-i <inputfile>] [-o <outputfile>] [-f <float>] [-p] [-s <imagefile>] [-q] [-y]')
      print('\t-p  display plot')
      print('\t-f  first (initial) value (default=0.0)')
      print('\t-s  save plot as imagefile (eps, jpeg, jpg, pdf, pgf, png, ps, raw, rgba, svg, svgz, tif, tiff)')
      print('\t-i  inputfile (default=/dev/stdin)')
      print('\t-o  outputfile (default=/dev/stdout)')
      print('\t-q  quiet; suppress output')
      print('\t-y  include y-data in output')
      sys.exit()
    elif opt in ("-i", "--ifile"):
      ifile=arg
      fin=open(arg,'r')
    elif opt in ("-o", "--ofile"):
      ofile=arg
    elif opt in ("-q"):
      ofile="/dev/null"
    elif opt in ("-d"):
      dmethod=int(arg)
    elif opt in ("-f"):
      f=float(arg)
    elif opt in ("-w"):
      w=int(arg)
    elif opt in ("-p"):
      plot=True
    elif opt in ("-s"):
      graph=True
      imagefile=arg
    elif opt in ("-y"):
      ydata=True

  data=np.loadtxt(ifile)
  x=data[:,0]
  y=data[:,1]

  y_int=integrate.cumtrapz(y,x,initial=0)

# Collect the data and write it
  if ydata:
      outdata=np.array([data[:,0],data[:,1],y_int])
      fmt=['%4.1f','%12.5e','%12.5e']
  else:
      outdata=np.array([data[:,0],y_int])
      fmt=['%4.1f','%12.5e']
  outdata=outdata.T
  np.savetxt(ofile,outdata,fmt)

# Build the plot

  fig=plt.figure()
  fig.set_size_inches(16,9)
  plt.plot(data[:,0],data[:,1],'-',color='black',label=ifile)
  plt.plot(data[:,0],y_int,'-',color='red',label=ofile)
  plt.legend()
  plt.minorticks_on()
  plt.grid(which='major',linestyle='-',linewidth='0.3')
  plt.grid(which='minor',linestyle=':',linewidth='0.3')
  if plot:
    plt.show()
  if graph:
    fig.savefig(imagefile,dpi=120)

main(sys.argv[1:])

我得到这个结果：

这里发生了什么？你知道吗

我将振幅为1的正弦波微分，得到振幅很小的余弦。当我把余弦和一个很小的振幅积分，得到一个振幅为1的正弦。余弦不也应该有振幅=1吗？你知道吗