Das deutsche Python-Forum

Code: Alles auswählen

from scipy import *
from scipy.optimize import leastsq
import csv
import Gnuplot

b = 0.00056
h = 0.00041
d = 0.000035

def residuals(p, y, x): 
	err = y-peval(x,p)
	return err

def peval(x, p):
	mm1 = []
	mm2 = []
	mm3 = []

	n = size(x,0)

	for i in range (0,n):
		z1 = min(p[0]/x[i], h/2)
		z2 = min(p[1]/x[i], h/2)
		
		m1 = 2*b*p[2]*(pow(z1,3)/3)*x[i]
		m2 = 2*b*(p[0]*(p[2]-p[3])*(pow(z2,2)-pow(z1,2))/2 + x[i]*p[3]*(pow(z2,3)-pow(z1,3))/3)
		m3 = 2.0*b* ((p[0]*(p[2]-p[3])+p[1]*(p[3]-p[4]))*(pow(h/2,2)-pow(z2,2))/2 + p[4]*(pow(h/2,3)-pow(z2,3))*x[i]/3)

		mm1.append(m1)
		mm2.append(m2)
		mm3.append(m3) 

# end for			

	M1 = array(mm1)
	M2 = array(mm2)
	M3 = array(mm3)	

	return M1+M2

e1_0 = 0.02
e2_0 = 0.07
E1_0 = 20000000000.0
E2_0 = 50000.0
E3_0 = 80000000000.0

x = []
y = []

file = ('Test.dat')
reader = csv.reader(open(file,'rb'),delimiter=';',skipinitialspace = True)
for row in reader:
	x.append(double(row[7])*1000.0) 
	y.append(double(row[6])/1000.0) 

# end for 

xx = array(x)
yy = array(y)

pname = (['e1','e2','E1','E2','E3'])
p0 = array([e1_0, e2_0, E1_0, E2_0, E3_0])
g = Gnuplot.Gnuplot()
d = Gnuplot.Data(x,y, title='Orig', with_='points 3 3')
g.xlabel('Curvature [1/mm]')
g.ylabel('Moment [Nmm]]')
g('unset grid')
plsq = leastsq(residuals, p0, args=(yy, xx), maxfev=10)
f = Gnuplot.Data(x,peval(x,plsq[0]),title='fit',with_='lines 1')
g.plot(f,d)


print "Final parameters"
for i in range(len(pname)):
	print "%s = %.4f " % (pname[i], p0[i])