Solving system of linear equations with different orders of magnitude

Question

I am trying to solve a system of two non-linear equations to obtain the value of one of the unknowns. The system of equations is:

p_sis=K*e^((-T/3)/(R_tot * C_tot))
p_dia=K*e^(-T/(R_tot * C_tot))

Data: p_sis = 22931.447 ; p_dia = 11865.691 ; T = 1.15384 ; R_tot = 9785025352

I have solved the system of equations through Wolfram Alpha and by hand and the results are the same:

K = 31865
C_tot = 1.18974*10^-10

I only need the value of C_tot. I would like to solve this system of equations automatically through Python. However, after having applied several possible methods, none of them have provided the correct answer (even when changing the initial guess values and using values very close to the actual solution).

Here are the different sets of code I have used:

METHOD #1: through scipy.optimize.fsolve

import numpy as np
from scipy.optimize import fsolve # 
p_sis = 22931.447
p_dia = 11865.691
T = 1.15384
R_tot = 9785025352 
def myFunction(z):
   x = z[0]  # = K
   y = z[1]  # = C_tot

   F = np.empty((2))
   F[0] = p_sis-x * np.exp((-T/3)/(R_tot * y))
   F[1] = p_dia-x * np.exp((-T)/(R_tot * y))
   return F

zGuess = np.array([30000,0.0000001])
z = fsolve(myFunction,zGuess)
print(z)

RESULT: [3.18787051e+04 4.14184981e-08].

METHOD #2: through simpy

from sympy import *
import numpy as np
import mpmath
p_sis = 22931.447
p_dia = 11865.691
T = 1.15384
R_tot = 9785025352 
mpmath.mp.dps=40
x = Symbol('x')  # = K
y = Symbol('y')  # = C_tot
aaa= nsolve([p_sis - x * 2.71828**(((-T / 3) / (R_tot * y))), p_dia - x * 2.71828**(((-T) / (R_tot * y)))], [x, y], [30000, 0.000000001])

print(aaa)

RESULT: Could not find root within the given tolerance.

METHOD #3: Least Squares Method with Levenberg-Marquardt (scipy.optimize.root)

import numpy as np
import sympy
import pandas as pd
from scipy.optimize import fsolve, root
def func(x):
    return [p_sis - x[0] * np.exp(((-T / 3) / (R_tot * x[1]))), p_dia - x[0] * np.exp(((-T) / (R_tot * x[1])))]

def jac(x):
    return np.array([[-np.exp(-T/(3*R_tot*x[1])),
                      -(T*x[0]*np.exp(-T/(3*R_tot*x[1])))/(3*R_tot*x[1]**2)],
                     [-np.exp(-T/(R_tot*x[1])),
                      -(T*x[0]*np.exp(-T/(R_tot*x[1])))/(R_tot*x[1]**2)]])


initialGuess = np.array([30000, 1e-10])
results = root(func, initialGuess, method='lm', jac=jac, tol=(1e-90), options={'maxiter': 5000000})

print(results)

RESULT: [2.89750393e+07, 8.07125772e-10].

METHOD #4: Newton's Method

from sympy import *
import numpy as np
p_sis = 22931.447
p_dia = 11865.691
T = 1.15384
R_tot = 9785025352 
def Newton_system(F, J, x, eps):
    from numpy import pi, exp
    F_value = F(x)
    F_norm = np.linalg.norm(F_value, ord=2)  # l2 norm of vector
    iteration_counter = 0
    while abs(F_norm) > eps and iteration_counter < 100:
        delta = np.linalg.solve(J(x), -F_value)
        x = x + delta
        F_value = F(x)
        F_norm = np.linalg.norm(F_value, ord=2)
        iteration_counter += 1
    # Here, either a solution is found, or too many iterations
    if abs(F_norm) > eps:
        iteration_counter = -1
    return x, iteration_counter

def test_Newton_system1():
    from numpy import pi, exp
    def F(x):
        return np.array(
            [p_sis - x[0] * exp(((-T / 3) / (R_tot * x[1]))),
             p_dia - x[0] * exp(((-T) / (R_tot * x[1])))])
    def J(x):
        return np.array(
            [[-exp(-T/(3*R_tot*x[1])), -(T*x[0]*exp(((-T / 3) / (R_tot * x[1]))))/(3*R_tot*x[1]**2)],
             [-exp(-T/(R_tot*x[1])), -(T*x[0]*exp(((-T) / (R_tot * x[1]))))/(R_tot*x[1]**2)]])
    expected = np.array([10000, 0.000001])
    tol = 1e-16
    x, n = Newton_system(F, J, x=np.array([0.02, 0.2]), eps=0.00000001)
    print(n, x)
    error_norm = np.linalg.norm(expected - x, ord=2)
    assert error_norm < tol, 'norm of error =%g' % error_norm
    print( 'norm of error =%g' % error_norm)
aaaaa=test_Newton_system1()
print(aaaaa)

RESULT: Singular matrix (did not reach a solution).

METHOD #5: Eq from sympy

from sympy import symbols, Eq, solve
import numpy as np
p_sis = 22931.447
p_dia = 11865.691
T = 1.15384
R_tot = 9785025352 
x, y = symbols("x y") # x = K ; y = C_tot ; e = approx 2.71828
equation_1 = Eq((p_sis - x * 2.71828**((-T/3)/(R_tot * y))), 0)
equation_2 = Eq((p_dia - x * 2.71828**((-T)/(R_tot * y))), 0)
print("Equation 1:", equation_1)
print("Equation 2:", equation_2)
solution = solve((equation_1, equation_2), (x, y))
print("Solution:", solution)

RESULT: computes forever and does not return a solution.

I have even conjoined the two equations into one and tried to solve it to obtain only C_tot:

p_sis*e^((T/3)/(R_tot * C_tot)) = p_dia*e^((T)/(R_tot * C_tot))

The code looks like this:

from math import exp
import scipy.optimize
p_sis = 22931.447
p_dia = 11865.691
T = 1.15384
R_tot = 9785025352 
    def fun(y):  # y = C_tot
    z = p_dia*exp((T)/(R_tot * y)) - p_sis*exp((T/3)/(R_tot * y))
    return z
z = scipy.optimize.fsolve(fun,10e-6, maxfev=1000)
print (z)

RESULT: The iteration is not making good progress, as measured by the improvement from the last ten iterations (solution equal to initial guess of 10e-6).

I suspect that the fact that K and C_tot have very different magnitudes, the program has trouble converging to the solution. Increasing the number of decimal places didn't really impact the performance.

How could I solve this issue? The fact that all the used methods did not reach the final solution makes me wonder if I made any mistakes in the definition of the problem.

Thank you in advance for all the help.