Python Programming
 
Forums: » Register « |  User CP |  Games |  Calendar |  Members |  FAQs |  Sitemap |  Support | 
User Name:
Password:
Remember me

The Shed is going Social! Join us on FaceBook and Twitter and chime in on the conversation.

Go Back   Dev Shed Forums > Programming Languages > Python Programming
Reload this Page

Problem with scipy.optimize.fsolve


Reply
Add This Thread To:
  Del.icio.us   Digg   Google   Spurl   Blink   Furl   Simpy   Y! MyWeb 
« Previous Thread | Next Thread »  
Thread Tools Search this Thread Rate Thread Display Modes
 
Unread Dev Shed Forums Sponsor:
  #1  
Old November 26th, 2012, 07:31 AM
leonardo2887 leonardo2887 is offline
Registered User
Dev Shed Newbie (0 - 499 posts)
  
Join Date: Nov 2012
Posts: 2 leonardo2887 User rank is Just a Lowly Private (1 - 20 Reputation Level) 
Time spent in forums: 44 m 41 sec
Reputation Power: 0
[Solved] Problem with scipy.optimize.fsolve
Dear all,

I am new to this forum as well as new to python. I started writing a program which I wanted to use to fit some experimental data with a particular function. In order to write the function, a non linear equation has to be solved.
In the few code lines below I am trying just to print out the function together with the independent variable (Ctot). I will then change the print command with a curve_fit() commnad.

However, to the program:

python Code:
Original - python Code 
  1.  
  2. import scipy
  3. import scipy.optimize, scipy.special, scipy.stats
  4. from math import sqrt
  5. from numpy import loadtxt
  6. from pylab import *
  7. from scipy.optimize import fsolve
  8.  
  9.  
  10.  
  11.        
  12. def Sig(Ctot, K, u, DHNC, DHC):
  13.     CP = 0.5
  14. #Determines the free surfactant concentration from the binding isotherm and mass balance
  15.     def func(x, Ctot, CP, K, u):
  16.          S = x-0.5*(1+(K*u*(Ctot-CP*x)-1)/sqrt((1-K*u*(Ctot-CP*x))**2+4*K*(Ctot-CP*x)))
  17.          #print "S=", S
  18.          return S
  19.     Theta = fsolve(func, 0.5, args=(Ctot, CP, K, u))
  20.     print "Theta =", Theta
  21.     Cf = Ctot-CP*Theta
  22.     lamb = 1. + sqrt((K*Cf*Theta)/(1.-Theta)) #lambda0, eigenvalue from weight matrix
  23.     chi = (Cf*(1.-Theta)*K)/(Theta*lamb**2) #fraction of non-cooperatively bound surfactant
  24.     s = K*u*Cf #new variable    
  25.     Theta1cf = 0.5*((u*K)/sqrt((s-1)**2+4*s/u)-((s-1)*(2*u*K*(s-1)+4*K))/(2*sqrt((s-1)**2+4*s/u)**3)) #first derivative of theta with respect to cf
  26.     a1 = -(4.0*Cf*Theta1cf*K)/(Theta*(sqrt((1-s)**2+4*s/u)+s+1)**2) #set of dummy functions to calculate S(Cf)
  27.     a2 = -(4.0*Cf*(1-Theta)*Theta1cf*K)/(Theta**2*(sqrt((1-s)**2+4*s/u)+s+1)**2)
  28.     a3 = +(4.0*(1-Theta)*K)/(Theta*(sqrt((1-s)**2+4*s/u)+s+1)**2)
  29.     a4 = -(8.0*Cf*(1-Theta)*K*((4*K-2*u*K*(1-s))/(2*sqrt((1-s)**2+4*s/u))+u*K))/(Theta*(sqrt((1-s)**2+4*s/u)+s+1)**3)
  30.     chi1cf = a1 + a2 + a3 + a4 #first derivative of chi with respect to Cf
  31.     DHm=DHNC*chi+(1-chi)*DHC #mean enthalpy of adsorption
  32.     return CP*(Theta1cf*DHm + Theta*(DHNC-DHC)*chi1cf)/(CP*Theta1cf+1)*0+Theta #calorimetric signal
  33.  #   return Theta
  34.     
  35. K = 2.
  36. u = 10.
  37. DHNC = -1000
  38. DHC = -3500
  39.  
  40. Ctot, Y = loadtxt('JR400-SDS-05.dat', unpack=True)
  41.  
  42.  
  43. print Ctot, Sig(Ctot, K, u, DHNC, DHC)
  44.  


The equation which has to be solve is contained in func(). Exactly the same output is obntained if func()... is put outside the function Sig(). If I run the program I get the following output:
Code:
[ 0.0203  0.028   0.0318  0.0355  0.0392  0.0429  0.0465  0.05    0.057
  0.0604  0.0671  0.0775  0.1155  0.1532  0.1905  0.2273  0.2638  0.2999
  0.3356  0.3709  0.4058  0.4404  0.4746  0.5084  0.5419  0.5751  0.6079
  0.6404  0.6726  0.7044  0.7359  0.7671  0.798   0.8286  0.8589  0.8889
  0.9186  0.948   0.9771  1.006   1.0346  1.0629  1.0909  1.1187  1.1462
  1.2431  1.6506] Theta = [ 0.0233803]
[ 0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803]


So, going with order, Ctot should be a ndarray. Also the output of Sig() is an ndarray with the same length as Ctot. What I really don't understand is, why is Theta an array of length 1?

I thought the problem may be, that Ctot is passed as an argument of the function, So I decided to print the values of the function (and uncommented the line print "S=", S in the function func() and the following output is optained:

Code:
[ 0.0203  0.028   0.0318  0.0355  0.0392  0.0429  0.0465  0.05    0.057
  0.0604  0.0671  0.0775  0.1155  0.1532  0.1905  0.2273  0.2638  0.2999
  0.3356  0.3709  0.4058  0.4404  0.4746  0.5084  0.5419  0.5751  0.6079
  0.6404  0.6726  0.7044  0.7359  0.7671  0.798   0.8286  0.8589  0.8889
  0.9186  0.948   0.9771  1.006   1.0346  1.0629  1.0909  1.1187  1.1462
  1.2431  1.6506] S= [ 0.5153606   0.51571418  0.51589425  0.51607325  0.51625595  0.51644247
  0.51662769  0.51681139  0.51718988  0.5173792   0.51776316  0.51838895
  0.52100977  0.52409536  0.52682602  0.52296529  0.45442974  0.00158271
 -0.32610579 -0.41085581 -0.44222299 -0.45773801 -0.46679046 -0.4726678
 -0.47678152 -0.47981323 -0.48212972 -0.48396008 -0.4854415  -0.48666077
 -0.48768421 -0.48855512 -0.489305   -0.48995722 -0.49052955 -0.4910357
 -0.49148641 -0.49189024 -0.49225406 -0.49258457 -0.4928851  -0.4931595
 -0.49341099 -0.49364309 -0.49385718 -0.49450905 -0.49620455]
S= [ 0.5153606   0.51571418  0.51589425  0.51607325  0.51625595  0.51644247
  0.51662769  0.51681139  0.51718988  0.5173792   0.51776316  0.51838895
  0.52100977  0.52409536  0.52682602  0.52296529  0.45442974  0.00158271
 -0.32610579 -0.41085581 -0.44222299 -0.45773801 -0.46679046 -0.4726678
 -0.47678152 -0.47981323 -0.48212972 -0.48396008 -0.4854415  -0.48666077
 -0.48768421 -0.48855512 -0.489305   -0.48995722 -0.49052955 -0.4910357
 -0.49148641 -0.49189024 -0.49225406 -0.49258457 -0.4928851  -0.4931595
 -0.49341099 -0.49364309 -0.49385718 -0.49450905 -0.49620455]
S= [ 0.5153606   0.51571418  0.51589425  0.51607325  0.51625595  0.51644247
  0.51662769  0.51681139  0.51718988  0.5173792   0.51776316  0.51838895
  0.52100977  0.52409536  0.52682602  0.52296529  0.45442974  0.00158271
 -0.32610579 -0.41085581 -0.44222299 -0.45773801 -0.46679046 -0.4726678
 -0.47678152 -0.47981323 -0.48212972 -0.48396008 -0.4854415  -0.48666077
 -0.48768421 -0.48855512 -0.489305   -0.48995722 -0.49052955 -0.4910357
 -0.49148641 -0.49189024 -0.49225406 -0.49258457 -0.4928851  -0.4931595
 -0.49341099 -0.49364309 -0.49385718 -0.49450905 -0.49620455]
S= [ 0.51536061  0.51571419  0.51589426  0.51607325  0.51625596  0.51644248
  0.5166277   0.5168114   0.51718988  0.51737921  0.51776317  0.51838895
  0.52100978  0.52409537  0.52682603  0.52296529  0.45442977  0.00158278
 -0.32610577 -0.4108558  -0.44222298 -0.457738   -0.46679045 -0.47266779
 -0.47678151 -0.47981322 -0.48212972 -0.48396007 -0.48544149 -0.48666076
 -0.4876842  -0.48855512 -0.48930499 -0.48995721 -0.49052954 -0.49103569
 -0.4914864  -0.49189023 -0.49225405 -0.49258456 -0.49288509 -0.49315949
 -0.49341099 -0.49364308 -0.49385717 -0.49450905 -0.49620454]
S= [-0.26435873 -0.38752518 -0.4518721  -0.51296691 -0.57023555 -0.62222561
 -0.66713794 -0.70543147 -0.76762254 -0.79190657 -0.83083645 -0.87403956
 -0.9478984  -0.97541116 -0.98909115 -0.9971163  -1.00237378 -1.00606038
 -1.0087817  -1.01086954 -1.01252017 -1.01386032 -1.01496646 -1.01589448
 -1.01668605 -1.017369   -1.0179624  -1.01848412 -1.01894631 -1.0193573
 -1.01972617 -1.02005902 -1.02036081 -1.02063567 -1.020887   -1.02111766
 -1.02133006 -1.02152626 -1.02170801 -1.02187739 -1.02203508 -1.02218222
 -1.02231981 -1.02244918 -1.0225706  -1.02295327 -1.02405022]
S= [ 0.17850246  0.17856586  0.1784605   0.17824108  0.17787755  0.17733401
  0.17658942  0.17561133  0.17264316  0.17054945  0.16462252  0.14785745
 -0.17734815 -0.63439506 -0.75038486 -0.7880122  -0.80539998 -0.81515722
 -0.82133763 -0.82558041 -0.82866368 -0.83100704 -0.8328412  -0.83431444
 -0.8355262  -0.83653991 -0.83739768 -0.83813474 -0.8387747  -0.83933378
 -0.83982772 -0.84026718 -0.84066063 -0.84101487 -0.84133542 -0.84162682
 -0.84189282 -0.84213657 -0.8423607  -0.84256816 -0.84276006 -0.84293807
 -0.84310361 -0.84325846 -0.8434031  -0.84385452 -0.8451126 ]
S= [ 0.10095452  0.09564024  0.09175446  0.08685868  0.08055503  0.07244711
  0.06234627  0.04987657  0.01436871 -0.009486   -0.072742   -0.21701886
 -0.6964262  -0.82038206 -0.85954789 -0.87728228 -0.88718851 -0.89343088
 -0.89770125 -0.90079686 -0.90313922 -0.90497578 -0.9064492  -0.90765662
 -0.90866632 -0.90952283 -0.91025624 -0.91089291 -0.91145066 -0.91194173
 -0.91237859 -0.91276967 -0.91312175 -0.91344032 -0.9137299  -0.91399422
 -0.91423643 -0.91445913 -0.91466456 -0.91485526 -0.91503215 -0.91519665
 -0.91534998 -0.91549373 -0.91562827 -0.91604991 -0.917239  ]
S= [-0.16762416 -0.27002632 -0.33030406 -0.3923348  -0.45484414 -0.51514617
 -0.56972298 -0.61780217 -0.69804371 -0.7297723  -0.78062969 -0.83639141
 -0.92744069 -0.95915022 -0.97433616 -0.98304749 -0.98867162 -0.9925751
 -0.99543487 -0.99761632 -0.99933314 -1.00072193 -1.00186475 -1.00282113
 -1.00363513 -1.00433613 -1.00494423 -1.00547813 -1.00595052 -1.00637012
 -1.00674634 -1.00708551 -1.0073928  -1.00767245 -1.00792799 -1.00816237
 -1.00837808 -1.00857722 -1.00876161 -1.00893338 -1.00909322 -1.0092423
 -1.00938167 -1.00951266 -1.00963557 -1.01002265 -1.01113017]
S= [  4.76482953e-02   2.97678644e-02   1.68644513e-02   7.33096042e-04
  -1.98074383e-02  -4.57380447e-02  -7.70691621e-02  -1.13990120e-01
  -2.07421266e-01  -2.60821382e-01  -3.73299783e-01  -5.35616150e-01
  -8.10004282e-01  -8.78532934e-01  -9.04728036e-01  -9.17938564e-01
  -9.25805375e-01  -9.30973544e-01  -9.34613166e-01  -9.37308281e-01
  -9.39380972e-01  -9.41026981e-01  -9.42361212e-01  -9.43463877e-01
  -9.44392526e-01  -9.45185056e-01  -9.45867199e-01  -9.46462027e-01
  -9.46985175e-01  -9.47447387e-01  -9.47859846e-01  -9.48230106e-01
  -9.48564267e-01  -9.48867310e-01  -9.49143340e-01  -9.49395772e-01
  -9.49627471e-01  -9.49840852e-01  -9.50037973e-01  -9.50221208e-01
  -9.50391382e-01  -9.50549816e-01  -9.50697661e-01  -9.50836403e-01
  -9.50966388e-01  -9.51374516e-01  -9.52531912e-01]
S= [-0.03356636 -0.08413044 -0.11903666 -0.16038587 -0.20910712 -0.26443916
 -0.3229645  -0.38206802 -0.49734822 -0.54827914 -0.63453103 -0.7317376
 -0.87962475 -0.9236086  -0.94283128 -0.95328665 -0.95981079 -0.96423383
 -0.96741938 -0.96981812 -0.97168695 -0.97318646 -0.97441219 -0.97543227
 -0.97629641 -0.97703759 -0.97767831 -0.97823911 -0.97873398 -0.97917249
 -0.97956482 -0.97991784 -0.98023712 -0.98052722 -0.98079192 -0.98103439
 -0.98125727 -0.9814628  -0.98165291 -0.98182983 -0.98199432 -0.98214762
 -0.9822908  -0.98242529 -0.98255139 -0.98294799 -0.98407812]
S= [ 0.00642834 -0.0268322  -0.05048608 -0.07947446 -0.1152534  -0.15839015
 -0.2073521  -0.26070825 -0.37731261 -0.43430441 -0.53793338 -0.66269849
 -0.85319019 -0.9055521  -0.92728945 -0.93877174 -0.94580673 -0.95051724
 -0.95387972 -0.9563948  -0.95834411 -0.95990172 -0.96117065 -0.9622237
 -0.96311367 -0.96387546 -0.96453282 -0.96510733 -0.9656136  -0.96606168
 -0.96646215 -0.96682215 -0.96714745 -0.9674428  -0.9677121  -0.96795861
 -0.96818507 -0.9683938  -0.96858676 -0.96876625 -0.96893306 -0.96908845
 -0.96923353 -0.96936975 -0.96949744 -0.96989874 -0.97104002]
S= [  7.48916438e-04  -3.48837053e-02  -6.01406021e-02  -9.09675784e-02
  -1.28792603e-01  -1.74032082e-01  -2.24868399e-01  -2.79634740e-01
  -3.97108288e-01  -4.53518931e-01  -5.54699510e-01  -6.74836917e-01
  -8.57629146e-01  -9.08503044e-01  -9.29805526e-01  -9.41113131e-01
  -9.48062147e-01  -9.52724627e-01  -9.56057706e-01  -9.58553531e-01
  -9.60489564e-01  -9.62037612e-01  -9.63299446e-01  -9.64347093e-01
  -9.65232834e-01  -9.65991251e-01  -9.66645898e-01  -9.67218170e-01
  -9.67722582e-01  -9.68169102e-01  -9.68568246e-01  -9.68927107e-01
  -9.69251432e-01  -9.69545928e-01  -9.69814482e-01  -9.70060336e-01
  -9.70286215e-01  -9.70494422e-01  -9.70686922e-01  -9.70865997e-01
  -9.71032426e-01  -9.71187476e-01  -9.71332251e-01  -9.71468192e-01
  -9.71595621e-01  -9.71996159e-01  -9.73135638e-01]
S= [ -1.90751672e-05  -3.59756293e-02  -6.14500778e-02  -9.25247598e-02
  -1.30622514e-01  -1.76137918e-01  -2.27214523e-01  -2.82154814e-01
  -3.99711351e-01  -4.56032066e-01  -5.56875738e-01  -6.76405949e-01
  -8.58208639e-01  -9.08890700e-01  -9.30136771e-01  -9.41421630e-01
  -9.48359424e-01  -9.53015625e-01  -9.56344855e-01  -9.58838157e-01
  -9.60772449e-01  -9.62319243e-01  -9.63580146e-01  -9.64627083e-01
  -9.65512268e-01  -9.66270243e-01  -9.66924532e-01  -9.67496511e-01
  -9.68000678e-01  -9.68446993e-01  -9.68845962e-01  -9.69204674e-01
  -9.69528871e-01  -9.69823254e-01  -9.70091710e-01  -9.70337478e-01
  -9.70563280e-01  -9.70771420e-01  -9.70963858e-01  -9.71142878e-01
  -9.71309258e-01  -9.71464263e-01  -9.71608998e-01  -9.71744902e-01
  -9.71872298e-01  -9.72272735e-01  -9.73411976e-01]
S= [  5.51917419e-08  -3.59484214e-02  -6.14174490e-02  -9.24859642e-02
  -1.30576937e-01  -1.76085493e-01  -2.27156151e-01  -2.82092157e-01
  -3.99646722e-01  -4.55969707e-01  -5.56821784e-01  -6.76367066e-01
  -8.58194262e-01  -9.08881076e-01  -9.30128545e-01  -9.41413968e-01
  -9.48352040e-01  -9.53008397e-01  -9.56337723e-01  -9.58831087e-01
  -9.60765423e-01  -9.62312248e-01  -9.63573174e-01  -9.64620129e-01
  -9.65505328e-01  -9.66263313e-01  -9.66917612e-01  -9.67489598e-01
  -9.67993770e-01  -9.68440091e-01  -9.68839065e-01  -9.69197780e-01
  -9.69521980e-01  -9.69816366e-01  -9.70084824e-01  -9.70330595e-01
  -9.70556399e-01  -9.70764540e-01  -9.70956980e-01  -9.71136001e-01
  -9.71302382e-01  -9.71457389e-01  -9.71602124e-01  -9.71738029e-01
  -9.71865426e-01  -9.72265865e-01  -9.73405113e-01]
S= [  4.05550246e-12  -3.59484999e-02  -6.14175431e-02  -9.24860761e-02
  -1.30577069e-01  -1.76085644e-01  -2.27156319e-01  -2.82092338e-01
  -3.99646908e-01  -4.55969887e-01  -5.56821940e-01  -6.76367178e-01
  -8.58194303e-01  -9.08881104e-01  -9.30128568e-01  -9.41413990e-01
  -9.48352062e-01  -9.53008418e-01  -9.56337743e-01  -9.58831108e-01
  -9.60765443e-01  -9.62312268e-01  -9.63573194e-01  -9.64620149e-01
  -9.65505348e-01  -9.66263333e-01  -9.66917632e-01  -9.67489618e-01
  -9.67993790e-01  -9.68440111e-01  -9.68839085e-01  -9.69197800e-01
  -9.69522000e-01  -9.69816386e-01  -9.70084844e-01  -9.70330615e-01
  -9.70556418e-01  -9.70764559e-01  -9.70957000e-01  -9.71136021e-01
  -9.71302402e-01  -9.71457409e-01  -9.71602144e-01  -9.71738049e-01
  -9.71865445e-01  -9.72265885e-01  -9.73405132e-01]
S= [ -3.46944695e-18  -3.59484999e-02  -6.14175431e-02  -9.24860761e-02
  -1.30577069e-01  -1.76085644e-01  -2.27156319e-01  -2.82092338e-01
  -3.99646908e-01  -4.55969887e-01  -5.56821940e-01  -6.76367178e-01
  -8.58194303e-01  -9.08881104e-01  -9.30128568e-01  -9.41413990e-01
  -9.48352062e-01  -9.53008418e-01  -9.56337743e-01  -9.58831108e-01
  -9.60765443e-01  -9.62312268e-01  -9.63573194e-01  -9.64620149e-01
  -9.65505348e-01  -9.66263333e-01  -9.66917632e-01  -9.67489618e-01
  -9.67993790e-01  -9.68440111e-01  -9.68839085e-01  -9.69197800e-01
  -9.69522000e-01  -9.69816386e-01  -9.70084844e-01  -9.70330615e-01
  -9.70556418e-01  -9.70764559e-01  -9.70957000e-01  -9.71136021e-01
  -9.71302402e-01  -9.71457409e-01  -9.71602144e-01  -9.71738049e-01
  -9.71865445e-01  -9.72265885e-01  -9.73405132e-01]
Theta = [ 0.0233803]
[ 0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803
  0.0233803  0.0233803  0.0233803  0.0233803  0.0233803]


So, also S, which is calculated in func() is an array of the same length of Ctot. What do I do wrong? Why can't I get Theta as an array with the same length of Ctot, with the value of the n-th element resulting from the n-th value of Ctot?

I am sorry if these questions are trivial, but I can't really understand what I am doing wrong.

With best regards,
Leonardo
Reply With Quote
  #2  
Old November 26th, 2012, 08:57 AM
leonardo2887 leonardo2887 is offline
Registered User
Dev Shed Newbie (0 - 499 posts)
  
Join Date: Nov 2012
Posts: 2 leonardo2887 User rank is Just a Lowly Private (1 - 20 Reputation Level) 
Time spent in forums: 44 m 41 sec
Reputation Power: 0
Problem solved with by substituting line 19

Theta = fsolve(func, 0.5, args=(Ctot, CP, K, u))

with

Theta = fsolve(func, Ctot/2., args=(Ctot, CP, K, u))

The size of theta is given by the size of the initial guess, and as I put only one number, the size of theta was only 1.

With best regards,
Leonardo
Comments on this post
b49P23TIvg agrees: Self solve. Yeah!!!
Reply With Quote
Reply

Viewing: Dev Shed Forums > Programming Languages > Python Programming > Problem with scipy.optimize.fsolve

« Previous Thread | Next Thread »

Developer Shed Advertisers and Affiliates



Thread Tools  Search this Thread 
Search this Thread:

Advanced Search
Display Modes  Rate This Thread 
Linear Mode Linear Mode
Rate This Thread:


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts
vB code is On
Smilies are On
[IMG] code is On
HTML code is Off
View Your Warnings | New Posts | Latest News | Latest Threads | Shoutbox
Forum Jump

Forums: » Register « |  User CP |  Games |  Calendar |  Members |  FAQs |  Sitemap |  Support | 
  
 


Powered by: vBulletin Version 3.0.5
Copyright ©2000 - 2013, Jelsoft Enterprises Ltd.

© 2003-2013 by Developer Shed. All rights reserved. DS Cluster - Follow our Sitemap