CODEWITHSUNDEEP
pythonmaxnumerical-methods
user35202
user35202

Reputation: 419

How to find multiple maxima of numerical time series

I am solving numerically (in python) a first order differential equation using Euler's method. I build the solution from time t=0 until some arbitrary time proceeding in time steps of size 0.05 say. One example of the resulting solution is in the figure below enter image description here

I would like to find the maxima that we see in that figure (as well as the times they occur) and store them in a dictionary. If there were just one maximum in the entire time range, I could use this code

y=[xinitial,viinitial,ainitial]            
t=0
maximum=-20000
maxai={}
h=0.05
ai=-2.1
for i in range(0,3701):
   dydt=computederivs(y)
   y = euler(y,dydt,h)
   t+=h
   if y[0]>maximum:
      maximum = y[0]
 maxai[ai]=maximum

Since I have multiple local maxima I have to detect them as I move in time t by somehow checking when the function goes down after climbing for a few time steps. I also need to store the maxima in a list which is the value to the dictionary key. I am imagining this is a common enough task that there must be well known ways to do it more or less simply?

Upvotes: 3

Views: 1524

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 25972

The event you should be looking for is when the first derivative v changes its sign. For a maximum then test if the second derivative a is negative.

 yp, y = y, euler(y,dydt,h)
 if yp[1]*y[1] < 0 and y[2] < 0:
      # do what you need to do to list the maximum

You could use linear interpolation to get a more precise ordinate for the maximum position.

If you want to get serious results, take a higher order method for the ODE integration.

Upvotes: 3

Lukasz Tracewski
Lukasz Tracewski

Reputation: 11377

I'd propose using find_peaks from scipy.signal module. This function takes a one-dimensional array and finds all local maxima by simple comparison of neighbouring values. Optionally, a subset of these peaks can be selected by specifying conditions for a peak’s properties.

Here's a code snippet to get you started:

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import find_peaks

Fs = 8000
f = 5
sample = 8000
x = np.arange(sample)
y = np.sin(2 * np.pi * f * x / Fs)
peaks = find_peaks(y)

plt.scatter(peaks[0], np.ones(f), c='red')
plt.plot(x, y)
plt.xlabel('sample(n)')
plt.ylabel('voltage(V)')
plt.show()

enter image description here

Your maxima:

print(peaks[0])
[ 400 2000 3600 5200 6800]

Upvotes: 5

Related Questions