[SciPy-User] scipy.optimize fmin error

Warren Weckesser warren.weckesser@enthought....
Tue Oct 6 13:58:42 CDT 2009


Also, it looks like you are missing some parentheses in your definition 
of the c array.  If you define c like this:

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

(note the parentheses around array([0.07, ..., 4.19])-0.07), you get the 
same answer as your matlab code: 0.6307

Warren


Warren Weckesser wrote:
> In your call to fmin, you have reversed the cr and tr args. It should be 
> args=(tr,cr) to match the definition of the residuals() function.
>
> Warren
>
> 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
>>
>> #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
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>   
>>     
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>   



More information about the SciPy-User mailing list