Villa Sega
Villa Sega

Reputation: 31

Least-squares polynomial fitting

How to fit this polynomial with leastsq? The aim is to get a graph optimized from the experimental and analytical values:

from scipy.integrate import quad
import pylab as py
import numpy as np

x_data is temperature, y_data is analytical and y1_data is experimental

x_data=[329.74403, 329.21733, 328.73927, 328.25111, 327.75969, 327.26852, 
326.7746,326.28142, 325.78471,325.28635,324.78976]
y_data=[]
y1_data=[1.55e-06, 1.82e-06, 1.93e-06, 1.17e-06, 1.93e-06, 1.79e-06, 1.31e- 
06, 1.75e-06,1.68e-06,1.39e-06, 1.69e-06]

Tavg = sum(x_data)/len(x_data)
B=0.33
m0 = 1.0
st = 0.1 * Tavg

def p(Tc):
    return 1.0/(np.sqrt(2.0*np.pi)*st)*np.exp(-((Tc- 
    Tavg)**2.0)/(2.0*st**2.0))


def F(Tc, Ti):
    if (Tc-Ti) <0:
       return 0
    else:
       return ((Tc-Ti)/Tc)**B*p(Tc)

this function gives out the analytical values

def M(t):
    for Ti in x_data:
       k = quad(F, 0.0, 2000.0, args=(Ti,))[0]
       y_data.append(k*m0)

#graph of temperature vs analytical is done here
def plt():
    py.xlabel('Temperature')
    py.ylabel('Magnetisation')
    py.plot(x_data,y_data)
    py.savefig('graf.png')

I need to find a graph optimized using leastsq which minimizes the error between the analytical and experimental value. I have looked online for leastsq examples but I had a hard time understanding and applying it to my code. Any help and insight is welcome.

Upvotes: 1

Views: 4523

Answers (1)

ForceBru
ForceBru

Reputation: 44838

You can use numpy.polyfit to do the fitting and numpy.polyval to get the data to plot.

coefficients = numpy.polyfit(x_data, y_data, degree)

fitted_data = numpy.polyval(coefficients, x_data)

Example usage

Generate and plot some random data that looks like stock price data:

from pylab import *

data = 10 + np.cumsum(np.random.normal(0.1, 0.1, size=100))
plot(data); grid(True); show()

You get something like this:

enter image description here

Then, do the fitting (get the coefficients of a polynomial that approximates your data) and the data to plot (evaluate the polynomial given by the coefficients you got):

X = np.arange(0, data.size)

coeff = np.polyfit(X, data, 5)
Y_fitted = np.polyval(coeff, X)

plot(Y_fitted); grid(True); show()

The result looks like this:

enter image description here

But to see the power of np.polyfit, we need one more graph: original data (orange) vs the polynomial (blue):

enter image description here

As you can see, on this concrete interval the polynomial approximates your data pretty well, and the higher the degree of the polynomial, the better the approximation.

Upvotes: 2

Related Questions