0

I am still new with python and I have a problem wit curve fitting. The following program is a simplification of a bigger program that I create but it represent the problem that I have.

The problem is that I have a function which I called burger that I cannot fit a curve. This line : y=np.sqrt(y) : is a problem. When I remove it, i can fit it perfectly but that not the function I want.

How Can I do a fitting of this function y=np.sqrt(y)?

# -*- coding: utf-8 -*-
"""
Created on Wed Dec 11 22:14:54 2013

@author: 
"""
import numpy as np
import matplotlib.pyplot as plt
import pdb
import scipy.optimize as optimization
from math import *
from scipy.optimize import curve_fit

import math
import moyenne

####################Function Burger###############################

def burger(t, E1, E2, N,tau):
    nu=0.4 #Coefficient de Poisson
    P=50 #Peak force
    alpha=70.3 #Tip angle
    y=((((pi/2.)*P*(1.-nu**2.))/(tan(alpha)))*(1./E1 + 1./E2*(1.-np.exp(-t/tau)) + 1./((N)*(1.-nu))*t))
    y=np.sqrt(y)
    return y

#######exemple d'utilisation##########   
xlist=np.linspace(0,1,100) 
ylist=[ burger(t,3, 2,1,0.1) for t in xlist] 

#pdb.set_trace()
pa,j = curve_fit(burger,xlist,ylist)

yfit=[burger(x,*pa) for x in xlist]

plt.figure()
plt.plot(xlist,ylist,marker='o')
plt.plot(xlist,yfit)
plt.show()

1 Answers1

1

So, this probably won't be the best answer you get, but while you wait for others here are some things to think about.

First, since you are new to python maybe you don't know, or maybe there is reason to solve these things in list comprehension, but I don't think you need the list comprehensions. You can use the numpy math operations to handle a whole array at a time. Instead of

    y=((((pi/2.)*P*(1.-nu**2.))/(tan(alpha)))* ...

You can write

    y = ((((np.pi/2.)*P*(1.-nu**2.))/(np.tan(alpha)))* ...

Then instead of

    [ burger(t, 3., 2., 1., 0.1) for t in xlist] 

you can do

    burger(xlist, 3., 2., 1., 0.1)

This is will be a lot faster when you are working with arrays.

Secondly, just looking through a couple of things that were happening in the algorithm. It wasn't looking for your parameters in the right ranges. I looked up the algorithm it is using on the scipy.optimize page (here) and wikipedia says that the convergence is dependent on the initial guess and also that it finds the local, not global, minima (Sometimes your code hit negative values for the parameters which made the sqrt of y undefined for some cases). If there is a way you can give it a good initial guess then it should work ([1., 3., 3., 2] worked for me). My command that solved it was: pa,j = curve_fit(burger,xlist,ylist, [1., 3., 3., 2], maxfev=10000)).

Thirdly, the first error I got when I used your code was that it reached the max number of fevals. Add maxfev=10000 (or more if you need) as the last argument to curve_fit.

Check it out. If you can give your bigger problem an initial guess then maybe you'll get it to converge. Otherwise maybe a different algorithm could be more suitable?

Update: See this question for a more detailed explanation of why this works, but you can get it to work without a guess if you give it another kwg, diag.

Use:

    pa,j = curve_fit(burger,xlist,ylist, diag=(1./xlist.mean(), 1./ylist.mean()), maxfev=10000)
Community
  • 1
  • 1
cc7768
  • 327
  • 4
  • 10