1

I'm trying to solve and reproduce a solution of a system of non-linear equations using scipy. I've learned a lot from different posts, and I could figure out some issues I had before. It seems, however, based on the solution I posted above to my own question, that scipy.optimize is quite sensitive to the different available methods I'm trying to use nonlinear solvers as someone else pointed this here but I'm getting the same solutions, which I suppose are wrong. I say this because in the scientific paper I'm following (and trying to reproduce the results) they solved the same system with Mathematica, and got a different result which is more close to what one would expect. In the past, I only posted a minimal work example of my code, but now I decided to post the real thing.

I was wondering if Mathematica could easily manage the choice of the solver and tackle the problem. I initially wanted to preferably use Python, but right now my priority is to solve this problem, no matter the language. If someone figures out a way to solve with scipy much better.

Follow my code:

import numpy as np
from scipy import optimize


c1 = 43.524e-3 

Z1 = 1 

Z2 = 3  

b = 1.7e-10 

q = 1.6e-19

k_B = 1.38e-23

T = 296.15

N = 6.02e23

lB=0.71e-9

xi = 4.2

v1 = 4*np.pi*np.exp(1.0)*N*(1.0+Z1)*(xi-Z1**(-1.0))*b**(3.0)
v2 = 4*np.pi*np.exp(1.0)*N*(1.0+Z2)*(xi-Z2**(-1.0))*b**(3.0) 
epsilon =  q**2/(xi*k_B*T*b)
I = 46.311e-3


kappa2 = (3.29e7)*c1**(0.5) * (1e2)


def f(p,*args):
    theta1,theta2 = p
    c2 = args[0]
    eq1 = theta1-(c1*v1/10**3)*np.exp(-1-2*Z1*xi*(1-Z1*theta1-Z2*theta2)*np.log(1-np.exp(-kappa2*b)))
    eq2 =  theta2-c2*((v2/(10**3*np.exp(1)))*((10**3*theta1*np.exp(1))/(c1*v1))**3)   
    return(eq1,eq2)



c2 =  np.arange(0,150e-6,1e-6)

THETA1 = []
THETA2 =[]
THETA = []
for i in c2:

    sol=optimize.root(f,np.array([0.12024367, 0.00062351]),args=(i),method='hybr')
    theta1,theta2 = sol.x[0],sol.x[1]
    sol.message
    theta=3*theta2+theta1
    THETA1.append(theta1)
    THETA2.append(theta2) 
    THETA.append(theta)
    print (theta1,theta2,theta)

What is supposedly wrong is that theta1 drops down extremelly fast, the numerical values for the two first interaction are: 

theta1,theta2,theta
0.12154778645003049, 6.444608966498463e-30, 0.12154778645003049
6.961522754431233e-05, 0.17477403029874458, 0.5243917061237781

And here is an evolution of theta1, theta2 and theta as a function of c2.

enter image description here

EDIT Results from the aforementioned paper:

enter image description here

ziulfer
  • 1,339
  • 5
  • 18
  • 30
  • Did you print `sol.success` after each call of `root()`? – Warren Weckesser Apr 23 '19 at 22:41
  • I just did and got `True` – ziulfer Apr 23 '19 at 22:43
  • I ran your code, and then double-checked that `f(sol.x, i)` was, in fact, approximately (0, 0) for each i. If those results don't agree with the paper, you should double-check that you formulated the problem in the same way as the paper. – Warren Weckesser Apr 23 '19 at 22:50
  • My main concern is if Mathematica could be a more robust way to solve this. I've tried different options for `method=` and it seems that `scipy` is extremely sensitive to it. I was wondering how would be the case for an equivalent solver in Mathematica. – ziulfer Apr 23 '19 at 23:00
  • Nothing that you have shown suggests that the scipy solver is not working. Have you double-checked your translation from Mathematica to Python? You have two versions of `eq2` in the code--are you using the correct version? Can you include the Mathematica code in the question? – Warren Weckesser Apr 23 '19 at 23:16
  • (Ignore the comment about two versions of `eq2`--I just noticed that you edited it out.) – Warren Weckesser Apr 23 '19 at 23:24
  • I haven't done anything with Mathematica. I coded from scratch in Python and based on several other posts I thought that Mathematica is perhaps more robust in the sense that, I suppose, it has fewer options. Also, the initial guess in`scipy` is somehow a problem, it seems. I will look for the posts that report it. I had two versions of `eq2` but the lastest I'm using is the one I posted. I'll edit my question to avoid confusion. Thanks! – ziulfer Apr 23 '19 at 23:25
  • I don't think anyone can help you without more concrete evidence that something is wrong. – Warren Weckesser Apr 23 '19 at 23:30
  • Ok, I added the graph from the paper. But still, I understand that that doesn't add a lot to the discussion. – ziulfer Apr 23 '19 at 23:33
  • My guess is that there is something wrong with your translation of the formulas from the paper into Python, maybe even just a small typo. – Warren Weckesser Apr 23 '19 at 23:37
  • Could you provide a link to the paper you used? I hope that should help – Bracula Apr 29 '19 at 13:51

0 Answers0