CODEWITHSUNDEEP
pythonmatlabscipymathematical-optimization
Eduard
Eduard

Reputation: 3

scipy.optimize.least_squares and matlab lsqnonlin difference

I would like to find minimum of equations given in these script. It is looking very messy(but deep undestanding of equation is not needed- I suppose). At the end of def is the expression to minimize: 

 vys1=-Qd-40*sqrt(5)*sqrt((ch+cm)*ep*kb*Na*T)*w1
 vys2=fi0-fib-Q0/cq 
 vys3=fib-fid+Qd/cq1 
 vysf= np.array([vys1,vys2,vys3])
 return vysf

I write this script in matlab using lsqnonlin to compare the results. Matlab results seems much accurate. Result are (fi0,fib,fid)

Python 
[-0.14833481 -0.04824387 -0.00942132] Sum(value) ~1e-3.  
Matlab 
[-0,13253   -0,03253      -0,02131 ] Sum(value)~1e-15

Note that script has a check for typos in equation(if they are identical in python and matlab) for [fi0,fib,fid]=[-0.120, -0.0750 ,-0.011] the result are the same [vys1,vys2,vys3]-

python [0.00069376  0.05500097 -0.06179421]
matlab [0.0006937598,0.05500096 -0.06179421]

Are there any options in least_squares to improve results? Thanks for any help(sorry for misunderstanding english )

Python

import scipy as sc
import numpy as np
from math import sinh
import matplotlib as plt
from numpy import exp, sqrt
from scipy.optimize import leastsq,least_squares

def q(par,ep,Na,kb,T,e,gamaal,gamasi,gamax,k1,k2,k3,k4,cq,cq1,ch,cm):
    fi0,fib,fid=np.array([par[0],par[1],par[2]])
    AlOH= gamaal*k1*exp(e*fi0/(T*kb))/(ch + k1*exp(e*fi0/(T*kb)))
    AlOH2= ch*gamaal/(ch + k1*exp(e*fi0/(T*kb)))
    SiO= gamasi*k2*exp(e*fi0/(T*kb))/(ch + k2*exp(e*fi0/(T*kb)))
    SiOH= ch*gamasi/(ch + k2*exp(e*fi0/(T*kb)))
    X= gamax*k3*k4*exp(e*fib/(T*kb))/(ch*k4 + cm*k3 + k3*k4*exp(e*fib/  (T*kb)))
    XH= ch*gamax*k4/(ch*k4 + cm*k3 + k3*k4*exp(e*fib/(T*kb)))
    Xm= cm*gamax*k3/(ch*k4 + cm*k3 + k3*k4*exp(e*fib/(T*kb)))
    Q0=e*(0.5*(AlOH2+SiOH-AlOH-SiO)-gamax)
    Qb=e*(XH+Xm)
    Qd=-Q0-Qb
    w1=sc.sinh(0.5*e*fid/kb/T)
    vys1=-Qd-40*sqrt(5)*sqrt((ch+cm)*ep*kb*Na*T)*w1
    vys2=fi0-fib-Q0/cq 
    vys3=fib-fid+Qd/cq1 
    vysf= np.array([vys1,vys2,vys3])
    return vysf

kb=1.38E-23;T=300;e=1.6e-19;Na=6.022e23;gamaal=1e16;gamasi=1e16
gamax=1e18;k1=1e-4;k2=1e5;k3=1e-4;k4=1e-4;cq=1.6;cq1=0.2
cm=1e-3;ep=80*8.8e-12
ch1=np.array([1e-3,1e-5,1e-7,1e-10])

# Check the equations, if they are same
x0=np.array([-0.120,    -0.0750 ,-0.011])
val=q(x0,ep,Na,kb,T,e,gamaal,gamasi,gamax,k1,k2,k3,k4,cq,cq1,ch1[0],cm)
print(val)
w1=least_squares(q,x0, args=(kb,ep,Na,T,e,gamaal,gamasi,gamax,k1,k2,k3,
                             k4,cq,cq1,ch1[0],cm))
print(w1['x'])

matlab

function[F1,poten,fval]=test()
kb=1.38E-23;T=300;e=1.6e-19;Na=6.022e23;gamaal=1e16;gamasi=1e16;gamax=1e18;
k1=1e-4;k2=1e5;k3=1e-4;k4=1e-4;cq=1.6;cq1=0.2;ch=[1e-3];cm=1e-3;ep=80*8.8e-   12;
% Test if equation are same 
x0=[-0.120, -0.0750 ,-0.011];
F1=rovnica(x0,ch) ;
[poten,fval]= lsqnonlin(@(c) rovnica(c,ch(1)),x0);
function[F]=rovnica(c,ch) 
fi0=c(1);
fib=c(2);
fid=c(3);
aloh=exp(1).^(e.*fi0.*kb.^(-1).*T.^(-1)).*gamaal.*k1.*(ch+exp(1).^(e.* ...
fi0.*kb.^(-1).*T.^(-1)).*k1).^(-1);
aloh2=ch.*gamaal.*(ch+exp(1).^(e.*fi0.*kb.^(-1).*T.^(-1)).*k1).^(-1);
sioh=ch.*gamasi.*(ch+exp(1).^(e.*fi0.*kb.^(-1).*T.^(-1)).*k2).^(-1);
sio=exp(1).^(e.*fi0.*kb.^(-1).*T.^(-1)).*gamasi.*k2.*(ch+exp(1).^(e.* ...
fi0.*kb.^(-1).*T.^(-1)).*k2).^(-1);
Xm=cm.*gamax.*k3.*(cm.*k3+ch.*k4+exp(1).^(e.*fib.*kb.^(-1).*T.^(-1)) ...
.*k3.*k4).^(-1);
 XH=ch.*gamax.*k4.*(cm.*k3+ch.*k4+exp(1).^(e.*fib.*kb.^(-1).*T.^(-1)) ...
.*k3.*k4).^(-1);
Q0=e*(0.5*(aloh2+sioh-aloh-sio)-gamax);
Qb=e*(XH+Xm);
Qd=-Q0-Qb;
F=[-Qd+(-40).*5.^(1/2).*((ch+cm).*ep.*kb.*Na.*T).^(1/2).*sinh((1/2).*e.* ...
fid.*kb.^(-1).*T.^(-1));...
 fi0-fib-Q0/cq;...
(fib-fid+Qd/cq1)];
 end
end

Upvotes: 0

Views: 3444

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114781

There is a mistake in this line:

w1=least_squares(q,x0, args=(kb,ep,Na,T,e,gamaal,gamasi,gamax,k1,k2,k3,
                             k4,cq,cq1,ch1[0],cm))

You have the argument kb in the wrong spot. The signature of q is:

def q(par,ep,Na,kb,T,e,gamaal,gamasi,gamax,k1,k2,k3,k4,cq,cq1,ch,cm):

The argument kb is between Na and T. If you fix the args argument in the least_squares call:

w1 = least_squares(q, x0, args=(ep, Na, kb, T, e, gamaal, gamasi, gamax,
                                k1, k2, k3, k4, cq, cq1, ch1[0], cm))

then the output of the Python script is

[ 0.00069376  0.05500097 -0.06179421]
[-0.13253313 -0.03253254 -0.02131043]

which agrees with the Matlab output.

Upvotes: 3

Related Questions