[SciPy-User] scipy.optimize fmin error

Jochen Schroeder cycomanic@gmail....
Tue Oct 6 17:00:40 CDT 2009


Hi Oz,

just some small comments. 

On 10/06/09 20:24, Oz Nahum wrote:
> 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

It is usually considered better style to keep the namespace of modules,
i.e. instead of 
from numpy import *
do:
import numpy as np

same for pylab.

> 
> #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)]

AFAIK generally it is not a good idea to test for equality of floating
point numbers. Instead you should test that the difference is smaller
than a number epsilon, i.e. find(abs(c-cmax)<epsilon). This avoids
problems due to the limited numerical precision. However in your case
it's a lot easier to use argmax. The code above would then look like
this.
c=c.transpose()
icmax= c.argmax()
cmax = c[icmax]
tmax = t[icmax]
cr=c/cmax
crmax=cr[icmax]
tr=t/int(tmax)
trmax=tr[icmax]


this also results in a small speed up (although probably negligible for
your relatively small arrays).

Cheers
Jochen

> 
> 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
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user


More information about the SciPy-User mailing list