Python - Nelder Mead vs fsolve: pH equilibria optimization problem leads to numerical instability in Nelder Mead?
I'm creating a pH equilibrium calculator - starting with initial compositions of certain acids / bases the calculator should come up with the equilibrated concentrations of each component based on physical constants. I've formulated this as a set of equations which need to be solved simultaneously. One key constraint is that H_plus*OH_minus must be 10^-14 , and I've included this in the equations. Code #1: Fsolve. This works, but is sensitive to initial conditions, and I can't put constraints (I need concentrations to be >=0). import numpy as np from scipy.optimize import fsolve def equations(variables, Ka, Kw, initial_concentrations): HA, A_minus, H_plus, OH_minus, ion_Na = variables HA_init, A_minus_init, H_plus_init, OH_minus_init, ion_Na_init = initial_concentrations # Define the equations eq1 = HA + A_minus - HA_init - A_minus_init # Mass balance on acetate eq2 = H_plus + ion_Na - OH_minus - A_minus # charge balance eq3 = H_plus *...