[SciPy-User] scipy.optimize fmin error

Oz Nahum nahumoz@gmail....
Tue Oct 6 13:24:36 CDT 2009


Hi Guys,
I'm trying to convert a matlab code to python, but I'm having too much
difficulties.
It's almost two days I'm struggling with the differences between the
languages. My latest success is actually using the function fmin
(before that I had success with leastsq).
But I can't can't a reasonable result comparing to matlab or my
knowledge (or early assumptions on my observations).

I attach here the code.

I don't understand why alpha (the dispersivity) I'm trying to find is
always calculated to ~24m where in matlab I get a result of ~0.6m.
Note also the difference in quality of models (the hand picked values,
are the ones I got from matlab)

Here's my code:

from pylab import *
from numpy import *
from scipy.optimize import leastsq, fmin

#injected mass in Kg
M = 0.01;
#distance between  the wells
r = 8.81; #[m]
## pumping rate
Q = 0.0061; #[m^3/sec]
##thickness of aquifer saturated with water
b=4.61; #[m];
##uncertainty of the measurments  (concentration measurments)
sigma_s = 0.01; # [m]

####
## define the measurments
####

t = array([1.0,300.0,600.0,	900.,	1200.,	1260.,	1320.,	1380, \
	1440,	1500,	1560,	1620,	1680,	1740,	1800,	1860,
	1920,	1980,	2040,	2100,	2160,	2220,	2280,	2340,\
	2400,	2460,	2520,	2580,	2640,	2700,	2760,	2820,\
	2880,	2940,	3000,	3060,	3120,	3180,	3240,	3300,\
	3360,	3420,	3480,	3540,	3600,	3660,	3720,	3780,\
	3840,	3900,	4200,	4500,	4800,	5100,	5400,	5700,\
	6000,	6300,	6600,	6900,	7200,	7500,	7800,	8100,\
	8400,	8700,	9000,	9300,	9600,	9900,   10200,	10500,\
	10800,	11100,	11400,	12000])

t = t.transpose()

c = array([0.07,	0.1,	0.11,	0.13,	1.17,	2.15,	3.65,	5.64,\
	8.12,	11,	14.3,	17.3,	20.6,	23.5,	26.5,	29.1,\
	31.5,	33.5,	35.3,	36.8,	37.9,	38.8,	39.5,	39.8,\
	40.1,	40.2,	40.1,	39.9,	39.5,	39,	38.5,	37.9,	\
	37.3, 36.5,	35.9,	35.1,	34.4,	33.5,	32.9,	32,	\
	31.2,	30.5,	29.9,	29,	28.2,	27.5,	26.8,	26.1,	\
	25.4,	24.7,	21.7,	19,	16.8,	14.8,	13.3,	12.1,	\
	11,	10.1,	9.4,	8.81,	8.15,	7.71,	7.3,	6.98,	\
	6.67,	6.36,	6.12,	5.92,	5.78,	5.58,	5.41,	5.15,	\
	4.77,	4.54,	4.37,	4.19])-0.07*1e-9*1350

c=c.transpose()
cmax = max(c)
## die index mit find(c==max(c)) kann man auch finden
tmax = t[find(c==cmax)]# 2460 ; #sec
cr=c/cmax # dimensionless concentration [-]
crmax = max(cr)
tr=t/int(tmax) # dimensionless time [-]
trmax = tr[find(cr==crmax)]

def residuals(alpha, tr, cr):
	#defintion for Radial Flow Field
	P = r/alpha#Peclet Number [-]
	#using tmax here causes an error of the optimization
	#it should be like the matlab version
	tmax=sqrt(1+P**(-2))-1/P;
	#see eq. 21 in Sauty, 1980, this creates a tmax(P), and P(alpha)
	K = tmax**0.5*exp((P/4/tmax)*(1-tmax)**2)
	print K
	A=-P/(4*tr)*(1-tr)**2#dimensionless
	f=K/tr**0.5*exp(A)#dimensionless
	#err=(f-cr)
	B=(f-cr)
	B=B**2
	N=len(B)
	B=sum(B)/N
	return B

def OneDmodel(alpha, r):
	P = r/alpha# #Peclet Number [-]
	tmax=sqrt(1+P**(-2))-1/P# %see eq. 21 in Sauty, 1980, this creates a
tmax(P), and P(alpha)
	K = sqrt(tmax)*exp(P/(4*tmax)*(1-tmax)**2)
	print K
	A=-P/(4*tr)*(1-tr)**2#; %dimensionless
	f=K/tr**0.5*exp(A) #%dimensionless
	return f

p0 = 3 #initial alpha value
#x = arange(0,6e-2,6e-2/30)
#alpha = leastsq(residuals, p0, args=(cr, tr))
alpha = fmin(residuals, 28, args=(cr,tr),maxiter=10000, maxfun=10000)
#print plsq[0]
#print alpha
#print plsq
print 'optimized dispersivity is ', alpha[0]
alpha=alpha[0]
oneDmodel = OneDmodel(alpha,r)
### Plot
handpicked= OneDmodel(0.6307,r)
plot(tr,cr, 'r+-')
plot(tr, oneDmodel, 'bo-')
plot(tr, handpicked, 'g--')

#print cr
#,x,y_meas,'o',x,y_true)
legend(['Real', 'Fit', 'Hand Picked'])
show()


Any kind of help will be more than appreciated !
Thanks

Oz Nahum
Graduate Student
Zentrum für Angewandte Geologie
Universität Tübingen

---

Imagine there's no countries
it isn't hard to do
Nothing to kill or die for
And no religion too
Imagine all the people
Living life in peace


More information about the SciPy-User mailing list