2

I wand to do integration on a dynamic system like to the ODE

x_ddot + d*x_dot + k*x = a*sin(omega*t)

with one modification: The external force a * sin(omega * t) must be replaced by another periodic signal which is based on some measurement data. That means: I need do curve fitting on my discrete measurement data and the resulting function must be periodic. I had two ideas how to solve the problem:

1) Use Fourier Transform (numpy.fft). But its difficult to generate a continuous function from discrete data using discrete fourier transform.

2) Use curve fitting, with a function like a1+a2*sin(omega*t)+a3*sin(2*omega*t)+a4*sin(.... where omega is 2*pi/(length of measurement data). Unfortunately this wasn't very successful as well. I've tried it with combination of sin and cos terms and went to very high order (...sin(10*omega*t)) which didn't make it better.

Plot

I would really appreciate some hints, please see my code with example-data bellow:

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math as math
import scipy.special as sp

# create data
f =    [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1.1,1.3,1.7,1.9,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2.1,2.3,2.9,4.1,4.3,4.4,4.5,4.4,4,3.6,3.2,2.8,2.4,2.0,1.6,1.3,1.2,1.1,1,1,1,1,1]
alpha   = np.linspace(0,len(f),len(f))

omega = 2*np.pi/len(f)

def func(alpha, a1, ac2, ac3, ac4, ac5, ac6, ac7):
    fit     = a1 + ac2*np.sin(omega*alpha) + ac3*np.sin(omega*2*alpha) + ac4*np.sin(omega*3*alpha) + ac5*np.sin(omega*4*alpha) +      ac6*np.sin(omega*5*alpha) + ac7*np.sin(omega*6*alpha)
    return fit

popt, pcov = curve_fit(func, alpha, f) 

print popt

y_fit= func(alpha,popt[0],popt[1],popt[2],popt[3],popt[4],popt[5],popt[6])

plt.plot(alpha,f,'bo',label='discrete data')
plt.plot(alpha,y_fit,'r',label='periodic fit')
plt.legend()
plt.show()
PyNoob
  • 43
  • 1
  • 6

1 Answers1

3

If you are content with a numerical function, a simple approach based on interpolation is the following:

Step 1: Define an interpolation function based on the data set, which determines one period of the function. (I used the simple interp1d from Scipy, more elaborate alternatives are also available in Scipy.)

from  scipy.interpolate import interp1d    
f_interp_one_period = interp1d(alpha, f, 'cubic')

Step 2: Define a periodic function based on the above (I used as period the value 55 based on the data)

def f_periodic(alpha):
        return f_interp_one_period(numpy.mod(alpha,55))

The following plot show f_periodic over the range -100 to 200, and the original data

enter image description here

Stelios
  • 5,271
  • 1
  • 18
  • 32
  • Thank you for your answer! Unfortunately I think a numerical function is not satisfying my problem. I need something like f(t)=k+a*sin(b*t)+c*sin(d*t)+....because i'm using pydstool – PyNoob Jul 11 '16 at 16:28
  • @PyNoob In that case you can compute the [Fourier series approximation](http://mathworld.wolfram.com/FourierSeries.html) of the numerical function, which requires determination of the Fourier coefficients. A related post is [this](http://stackoverflow.com/questions/4258106/how-to-calculate-a-fourier-series-in-numpy) – Stelios Jul 11 '16 at 17:00
  • @Stelios consider fitting data as function of angle, going from 0 to 2pi. Is it possible to somehow impose f(0)=f(2pi)? Thank you. – Alexander Cska Jun 02 '20 at 15:10