#1
  1. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Nov 2012
    Posts
    2
    Rep 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:
     
    import scipy
    import scipy.optimize, scipy.special, scipy.stats
    from math import sqrt
    from numpy import loadtxt
    from pylab import *
    from scipy.optimize import fsolve
     
     
     
     
    def Sig(Ctot, K, u, DHNC, DHC):
        CP = 0.5
    #Determines the free surfactant concentration from the binding isotherm and mass balance
        def func(x, Ctot, CP, K, u):
             S = x-0.5*(1+(K*u*(Ctot-CP*x)-1)/sqrt((1-K*u*(Ctot-CP*x))**2+4*K*(Ctot-CP*x)))
             #print "S=", S
             return S
        Theta = fsolve(func, 0.5, args=(Ctot, CP, K, u))
        print "Theta =", Theta
        Cf = Ctot-CP*Theta
        lamb = 1. + sqrt((K*Cf*Theta)/(1.-Theta)) #lambda0, eigenvalue from weight matrix
        chi = (Cf*(1.-Theta)*K)/(Theta*lamb**2) #fraction of non-cooperatively bound surfactant
        s = K*u*Cf #new variable    
        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
        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)
        a2 = -(4.0*Cf*(1-Theta)*Theta1cf*K)/(Theta**2*(sqrt((1-s)**2+4*s/u)+s+1)**2)
        a3 = +(4.0*(1-Theta)*K)/(Theta*(sqrt((1-s)**2+4*s/u)+s+1)**2)
        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)
        chi1cf = a1 + a2 + a3 + a4 #first derivative of chi with respect to Cf
        DHm=DHNC*chi+(1-chi)*DHC #mean enthalpy of adsorption
        return CP*(Theta1cf*DHm + Theta*(DHNC-DHC)*chi1cf)/(CP*Theta1cf+1)*0+Theta #calorimetric signal
     #   return Theta
     
    K = 2.
    u = 10.
    DHNC = -1000
    DHC = -3500
     
    Ctot, Y = loadtxt('JR400-SDS-05.dat', unpack=True)
     
     
    print Ctot, Sig(Ctot, K, u, DHNC, DHC)


    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
  2. #2
  3. No Profile Picture
    Registered User
    Devshed Newbie (0 - 499 posts)

    Join Date
    Nov 2012
    Posts
    2
    Rep 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!!!

IMN logo majestic logo threadwatch logo seochat tools logo