[SciPy-User] scipy.optimize fmin error
Warren Weckesser
warren.weckesser at enthought.com
Tue Oct 6 14:58:42 EDT 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 at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-User
mailing list