scipy.optimize.least_squares и matlab lsqnonlin разница

Question

Я хотел бы найти минимум уравнений, приведенных в этом скрипте. Это выглядит очень грязно (но глубокая неопределенность уравнения не нужна - я полагаю). В конце def это выражение для минимизации:

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 

Я пишу этот скрипт в MATLAB с помощью lsqnonlin для сравнения результатов. Результаты Matlab кажутся очень точными. Результат является (fi0, выдумка, FID)

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

Обратите внимания, что скрипт имеет чек на опечатки в уравнении (если они идентичны в питоне и MATLAB) для [fi0,fib,fid]=[-0.120, -0.0750 ,-0.011] результата те же [vys1,vys2,vys3] -

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

Есть еще варианты в least_squares для улучшения результатов? Спасибо за любую помощь (извините за недопонимание английский)

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 

Answer 1

Существует ошибка в этой строке:

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

У вас есть аргумент kb в неправильном месте. Подпись q является:

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

Аргумент kb между Na и T. Если зафиксировать args аргумент в least_squares вызова:

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

, то выход из сценария Python является

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

, что согласуется с выходом Matlab.