1

I'm trying to minimise the function ellipmodel

import numpy as np
import matplotlib.pyplot as plt
import scipy as sp
import sympy as sym    

def ellipmodel(c,t,y,dt):
    m = (np.sin(np.exp(-c[1]*t)*c[0]/2))**2
    w = np.pi*c[2]/2/ellipke(m)
    phase = np.cumsum(w)*dt
    error = sum((y-c[0]*np.exp(-c[1]*t)*np.cos(phase+c[3])-c[4])**2)
    # The square of the difference between data and model is to be minimised.
    return error

in this function:

# c[0]=theta_m, c[1]=alpha, c[2]=omega_0,
# c[3]=phi, c[4]=const.
# pos is the position data

def expfit(c,pos):
    # Determines a fit to the large amplitude decay of the undriven pendulum.
    min = 150
    t = 0.001 * np.arange(min,len(pos)+1) # Converts data to seconds.
    pos = pos[min:len(pos)]
    ### options = optimset('MaxIter',5000) # Determines the number of iterations.
    dt = t[1] - t[0]
    # Minimise the function ellipmodel.
    ### c = fminsearch('ellipmodel',c(),options,t,pos,dt)
    m = (np.sin(np.exp(-c[1]*t)*c[0]/2))**2
    # Takes account of change of frequency with amplitude.
    w = np.pi*c[2]/2/sym.elliptic_k(m)
    phase = np.cumsum(w)*dt
    fit = c[0]*np.exp(-c[1]*t)*np.cos(phase+c[3])+c[4] # Uses c values from fmin.
    return fit

The fit parameters c are to be guessed, say c = np.array([1,1,1,1,1]). The pos array is data recorded in a .txt file (which I don't know how to upload), but you could try it with some function like pos = exp(-arange(0,6,0.01))*sin(0.5*arange(0,6,0.01)+3).

I'm actually translating this from Matlab to Python. I think I'm quuite close, but I'm not particularly sure about the lines above which start with a ###.

Ali
  • 56,466
  • 29
  • 168
  • 265
Medulla Oblongata
  • 3,771
  • 8
  • 36
  • 75
  • 1
    possible duplicate of [numpy/scipy analog of matlab's fminsearch](http://stackoverflow.com/questions/19070943/numpy-scipy-analog-of-matlabs-fminsearch) – horchler Feb 25 '14 at 20:22
  • @horchler: Yes, but is there a way to modify the answer given there to incorporate the `ellipmodel` function? – Medulla Oblongata Feb 25 '14 at 22:34
  • Of course there's a way. They're just different function names. You use `ellipmodel` and the other answer uses `banana`. I suggest reading the documentation for [`scipy.optimize.fmin`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin.html#scipy.optimize.fmin) and trying what the answers to the other question suggest. If you can't get it to work, then you might have a specific non-duplicate question about precisely what's not working. – horchler Feb 25 '14 at 22:51

0 Answers0