I'm moving from Excel to Python and I'm trying to solve these equations:

$$\begin{align} X_1&=\bigg[\big(3.47-\log(X_2)\big)^2+\big(\log(c)+1.22)^2\bigg]^{0.5}\\ X_2&=\frac{a}{101.32}\bigg(\frac{101.32}{b}\bigg)^{X_3}\\ X_3&=0.381X_1+0.05\bigg(\frac{b}{101.32}\bigg)-0.15 \end{align} $$ I have tried this code

import numpy as np
import pandas as pd
import scipy.optimize as opt
from scipy.optimize import NonlinearConstraint
a = np.linspace(300,400,30) 
b = np.random.randint(700,18000,30) 
c = np.random.uniform(1.4,4.0,30) 
df = pd.DataFrame({'A':a, 'B':b, 'C':c})

def func(zGuess,*Params):
    x,y,z = zGuess
    a,b,c = Params
    eq_1 = (((3.47-np.log10(y))**2+(np.log10(c)+1.22)**2)**0.5) - x
    eq_2 = ((a/101.32) * (101.32/b)** z) - y
    eq_3 = (0.381 * x + 0.05 * (b/101.32) -0.15) - z
    return eq_1,eq_2,eq_3

zGuess = np.array([2.6,20.2,0.92])
df['x'],df['y'],df['z'] = zip(*df.apply(lambda x: opt.fsolve(func,zGuess,args=(x['A'],x['B'],x['C'])),1) )

but it returns the first guess as solution.

In Excel I would run a VBA:

Sub Solver_n()
' Macro6 Macro

'Application.Calculation = xlCalculationManual
Application.ScreenUpdating = False
For l = 3 To 32

    SolverOk SetCell:=Cells(l, 14), MaxMinVal:=2, ValueOf:=0, ByChange:=Cells(l, 12), Engine:=1 _
        , EngineDesc:="GRG Nonlinear"
    SolverOk SetCell:=Cells(l, 14), MaxMinVal:=2, ValueOf:=0, ByChange:=Cells(l, 12), Engine:=1 _
        , EngineDesc:="GRG Nonlinear"
    SolverSolve (True)
Next l
'Application.Calculation = xlCalculationAutomatic
Application.ScreenUpdating = True
End Sub

where (l,14) is the residual and (l,12) is the n_calc and would give the correct solution as presented here.

I trust the excel solution because it has a physical meaning and it is what I expect, but the Python solution does not match. How can resolve this problem in Python?

I would like to move to Python because it is easier to deal with plotting and large datasets as Pandas dataframes.


1 Answer 1


The problem is that y gets values close to zero so the solver steps over the axis and tries negative values during the derivative computation via difference quotients. This does not work with log10.

You can

  • replace log10(y) with the cut-off version 0.5*log10(1e-8+y*y). This gives solution, you would have to check how good they are relative to the original equations.

  • replace y with 10**v, so that log10(y) is replaced with v and the second equation with

eq2 = (np.log10(a/101.32) -np.log10(101.32/b)*z) - v

You would have to translate the initial value of y with the one for v and the solution back,...

  • just iterate the given equations as fixed-point iteration. This might not converge, or not converge fast, but has no numerical singularities.

Some solvers also allow to set boundaries, so that negative values in y can be directly avoided.

  • $\begingroup$ Thank you! this answer helped me to understand better my problem! $\endgroup$
    – JCV
    Apr 4, 2022 at 9:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.