Python, Fitting into a system of equations

Question

I am trying to fit my data into a system of two equations like below:

Here I have gamma versus C_s, and rest of them, gamma_0, Gamma_inf, a, and K are fitting parameters.

I don't see any error, But the fitting does not work, and every parameters cannot be negative but some of them are negative.

import numpy as np
from scipy.optimize import curve_fit, fsolve
import matplotlib.pyplot as plt
import pandas as pd

# Prompt the user for the CSV file location
csv_file = input("Enter the path to your CSV file: ")

# Read the CSV file into a DataFrame
df = pd.read_csv(csv_file,header=None)

# Alternatively, select each column by index
# Remember that indices start from 0
Cs_data = df.iloc[:, 0]
gamma_data = df.iloc[:, 1]

# Define the Frumkin equilibrium function to solve for x
def frumkin_x(C_s, a, K, x_guess=0.5):
    # Function to find the root of
    func = lambda x: x - C_s / (C_s + a * np.exp(K * x))
    x_solution, = fsolve(func, x_guess)
    return x_solution

# Define the equilibrium isotherm function
def equilibrium_isotherm(C_s, gamma_0, Gamma_inf, a, K):
    x = np.array([frumkin_x(cs, a, K) for cs in C_s])
    return gamma_0 + Gamma_inf * R * T * (np.log(1 - x) - 0.5 * K * x**2)

# Constants
R = 8.314  # Universal gas constant, J/(mol*K)
T = 298    # Temperature, K

# Initial parameter guesses for curve_fit
# gamma_0, Gamma_inf, a, K
initial_guesses = [72, 0.0000000004, 0.00000007, 0]

params_opt, params_cov = curve_fit(lambda Cs, gamma_0, Gamma_inf, a, K: equilibrium_isotherm(Cs, gamma_0, Gamma_inf, a, K), Cs_data, gamma_data, p0=initial_guesses)

# Print optimized parameters
print("Optimized Parameters:")
print("gamma_0:", params_opt[0])
print("Gamma_inf:", params_opt[1])
print("a (alpha_0/beta_0):", params_opt[2])
print("K:", params_opt[3])

# Plot the original data and fitted curve for visualization
Cs_fit = np.linspace(min(Cs_data), max(Cs_data), 100)
gamma_fit = equilibrium_isotherm(Cs_fit, *params_opt)

plt.scatter(Cs_data, gamma_data, label='Data')
plt.plot(Cs_fit, gamma_fit, label='Fitted Curve', color='red')
plt.xlabel('Sublayer Concentration (Cs)')
plt.ylabel('Surface Tension (gamma)')
plt.legend()
plt.show()

I did not see any errors, so I cannot figure out what is the problem.

above link contains my script and data.

Answer 1

With a little bit of help, essentially bounds and initial guess (fsolve and curve_fit), you can regress your parameters:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import optimize

data = pd.read_csv("frumkin.csv", header=None, names=["Cs", "gamma"])

def system(x, Cs, a, K):
    return x - Cs / (Cs + a * np.exp(K * x))

R = 8.314
T = 298.

def model(Cs, g0, Ginf, a, K):
    x = np.array([optimize.fsolve(system, x0=0.05, args=(c, a, K))[0] for c in Cs])
    return g0 + Ginf * R * T * (np.log(1. - x) - 0.5 * K * x ** 2)

popt, pcov = optimize.curve_fit(
    model, data.Cs, data.gamma,
    p0=[100., 1e-5, 1e-5, 1.],
    bounds=[(0., 0., 0., 0.), (np.inf, np.inf, np.inf, np.inf)]
)
# array([7.28040792e+01, 3.58127986e-01, 1.43487197e-08, 3.72400451e+02])

It returns:

About uniqueness of the solution, one can show that it is a valid hypothesis over the parameter ranges. Following abacus show relationship between Cs and x for different setup:

@np.vectorize
def solution(Cs, a, K):
    return optimize.fsolve(system, x0=0.05, args=(Cs, a, K))[0]

Cs = np.logspace(-12, -6, 100, base=10)
As = np.logspace(-12, -6, 10, base=10)
Ks = np.logspace(-5, +5, 10, base=10)
C, A, K = np.meshgrid(Cs, As, Ks, indexing="ij")

s = solution(C, A, K)