Scipy Curve_fit: how would I go about improving this fit?

I've been working on a standard potential which I am trying to fit with a given model: ax² - bx³ + lx⁴

The x and y values for the fit are generated from the code as well, the x values are generated by numpy.linspace. I have bounded the a,b,c parameters such that they are always positive. I needed the fit to mimic the data in such a way that at minimum, the height of the local maxima and the position of the global minima are accurate. Instead this is what I'm getting (blue is the actual data and the dashed line is the fitted one with the given model):

Here is the relevant part of my code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.interpolate import *

v =246.
x = np.linspace(0.1, 270.1, 27001)  
Xdata = x/v

def func(x,a,b,l):
    return (a*(x**2)) - (b*(x**3)) + ((l)*(x**4))

temp = np.linspace(80.,80.,1)
b = np.zeros_like(temp)
a = np.zeros_like(temp)
l = np.zeros_like(temp)
Vfit = pot_giver(temp[0])   #External func.
Ydata = (Vfit - Vfit[0])/(pow(v,4))
popt, pcov = curve_fit(func, Xdata, Ydata,bounds=((0.,0.,0.), (np.inf,np.inf,np.inf)))
b[0],a[0],l[0] = popt

plt.plot(Xdata,Ydata)
plt.plot(Xdata, func(Xdata, *popt), 'g--')
plt.show()

I am doing this as an array because I am varying the temp over a wide range, this singular element 'temp' was created for troubleshooting.

The Xdata and Ydata are here on the first and second columns respectively: Data

You can also use scipy.optimize.leastsq

Here, you can also impose constraints (through a penalty added to the residuals when they are not satisfied).

You asked for the following constraints:

• the maximum of the function is correct;

• the position of the global minimum is correct.

It's a little sensitive to how you impose the constraints, but this just about gets those.

import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq

Xdata, Ydata = np.loadtxt( "Pot_data_online.txt", unpack=True )
ymax = Ydata.max()
ymin = Ydata.min()
xymin = Xdata[0]
for x, y in zip( Xdata, Ydata ):
    if y == ymin: xymin = x


# Fitting function
def f( x, a, b, l ):
    return a * x ** 2 - b * x ** 3 + l * x ** 4


# Residuals, including a penalty to get the constraints
def residuals( p, x, y ):
    a, b, l = p
    x1 = ( 3 * b - math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l )
    x2 = ( 3 * b + math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l )
    error = y - f( x, a, b, l )
    penalty = 1e6 * ( abs( f(x1,a,b,l) / ymax - 1 ) + abs( x2 / xymin - 1 ) )
    return abs( error ) + penalty

popt, pcov = leastsq( func=residuals, x0=(0.01,0.3,0.01), args=(Xdata,Ydata) )
a, b, l = popt
x1 = ( 3 * b - math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l )
x2 = ( 3 * b + math.sqrt( max( 9 * b ** 2 - 32 * a * l, 0 ) ) ) / ( 8 * l )
print( "a, b, l = ", a, b, l )
print( "Constraint on maximum:             ymax , fmax  = ", ymax , f( x1, a, b, l ) )
print( "Constraint on position of minimum: xymin, xfmin = ", xymin, x2               )
plt.plot( Xdata,Ydata, 'b-' )
plt.plot( Xdata, f( Xdata, *popt ), 'g--')
plt.show()

Output:

a, b, l =  0.026786218092281242 0.0754627163616853 0.04481075395203686
Constraint on maximum:             ymax , fmax  =  0.0007404005868381408 0.0007404009309691355
Constraint on position of minimum: xymin, xfmin =  0.947308 0.947621778120905