我正在研究一个使用二维傅立叶级数的优化问题,基本上我正在尝试实现以下公式:
我想在我的脚本中将其实现为阶数 m 的函数,在某种意义上说,只要提供 m,我就可以用 m 阶级数拟合一些数据,并获得级数的最佳系数(除了我已经知道的 C_00)
无论如何,当涉及到使用Scipy的curve_fit函数与参数数量可变的模型(例如,如果 m=1 我有8个系数+ 2个角速度,如果 m=2 我有24个系数+ 2个角速度等)时,我很挣扎。
到目前为止,我已经写了以下内容:
c00=50
def fourier(x, v, w1, w2):
f=c00
k=0
for i in range(1, m+1):
f=f+v[k]*np.cos(i*w1*x)+v[k+1]*np.sin(i*w1*x)
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+v[k]*np.cos(i*w1*x+j*w2*x)+v[k+1]*np.sin(i*w1*x+j*w2*x)
k+=2
return f
time = [] #some x data I get from a txt file
energy = [] # some y data I get from a txt file
popt, pcov = curve_fit(fourier, time, energy)
但当我运行脚本时,它引发了以下错误:
文件“/home/arcalinux/Desktop/fourier_series.py”,第63行,傅立叶格式f=f+v[k]np.cos(iw1x)+v[k+1]np.sin(iw1x)索引错误:标量变量的索引无效。
我还尝试了其他功能:
def fourier2(args):
f=c00
k=3
for i in range(1, m+1):
f=f+args[k]*np.cos(i*args[1]*args[0])+args[k+1]*np.sin(i*args[1]*args[0])
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+args[k]*np.cos(i*args[0]*x+j*args[1]*x)+
+args[k+1]*np.sin(i*args[0]*x+j*args[1]*x)
k*=2
return f
但是当我运行curve_fit时,我有:
ValueError:无法确定拟合参数的数目。
请注意,傅立叶级数函数正确工作,例如(随机数):
v=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
m=2
x=2
w1=3
w2=4
print(fourier(x,v,w1,w2))
49.136961626348295
同样的道理也适用于傅立叶2
我怎样才能解决这个问题呢?谢谢!
完整代码如下所示:
# importing libraries
import numpy as np
import math
from scipy.optimize import curve_fit
from statistics import mean
import matplotlib.pyplot as plt
m=2 #order of series
number_of_coeff=(2*m+1)**2 - 1
N=number_of_coeff
# build series
coefficients=[]
# first part, j=0
for i in range(1,m+1):
coefficients.append('c'+str(i)+'0')
coefficients.append('s'+str(i)+'0')
# second part
for j in range(1, m+1):
for i in range(-m, 0):
coefficients.append('g'+str(-i)+str(j))
coefficients.append('t'+str(-i)+str(j))
for i in range(0, m+1):
coefficients.append('c'+str(i)+str(j))
coefficients.append('s'+str(i)+str(j))
print(len(coefficients))
v=np.ones(N)
m=2
x=2
w1=3
w2=4
arg=[2,3,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
arg3=[3,4,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
c00=50
def fourier(x, v, w1, w2):
print(type(v))
print('v is ',v)
f=c00
k=0
print(v[0])
for i in range(1, m+1):
f=f+v[k]*np.cos(i*w1*x)+v[k+1]*np.sin(i*w1*x)
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+v[k]*np.cos(i*w1*x+j*w2*x)+
v[k+1]*np.sin(i*w1*x+j*w2*x)
k+=2
return f
def fourier2(args):
f=c00
k=3
for i in range(1, m+1):
f=f+args[k]*np.cos(i*args[1]*args[0])+
args[k+1]*np.sin(i*args[1]*args[0])
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+args[k]*np.cos(i*args[1]*args[0]+
j*args[2]*args[0])+
args[k+1]*np.sin(i*args[1]*args[0]+
j*args[2]*args[0])
k+=2
# print('number of args =',len(args))
return f
def fourier3(x,args):
f=c00
k=2
for i in range(1, m+1):
f=f+args[k]*np.cos(i*args[0]*x)+
args[k+1]*np.sin(i*args[0]*x)
k+=2
for j in range(1, m+1):
for i in range(-m, m+1):
f=f+args[k]*np.cos(i*args[0]*x+j*args[1]*x)+
args[k+1]*np.sin(i*args[0]*x+j*args[1]*x)
k+=2
# print('number of args =',len(args))
return f
# print(fourier2(arg))
# print(fourier3(x,arg3))
# print(fourier(2,v,3,4))
# getting data
time = []
energy = []
filepath=''
filename=''
with open(filepath+filename, 'r') as f:
for line in f: # iterate through the file
if not line: continue
t, e = line.split()
time.append(float(t)),
energy.append(float(e)) # add a value
c00=mean(energy)
time=np.array(time)
popt, pcov = curve_fit(fourier, time, energy)
print(popt)
1条答案
按热度按时间tyg4sfes1#
根据scipy documentation,您可以为
curve_fit
函数提供p0
参数,这是您对函数参数的初始猜测(这里是v
、w1
和w2
)。如果不这样做,函数将假定参数的初始值为1,这就是为什么
v
将得到标量而不是数组。因此,在调用
curve_fit
之前,您可以执行以下操作: