[FIXED] How do I make a function from titration points in python?

Issue

I want to automatically calculate a function by giving points on the graph, when I run the following script.

Titration curves are pretty similar to logistic regression, tho, they have different inclinations at the turning points. So it's not 100% the same curve. The best would be to have an output, that gives you the function, similar to this one: y = (A) / (1 + B^{x-C}) + D

Btw. sorry for my english, it's not my native language. Please ask questions, if it is not clear.

#from https://www.geeksforgeeks.org/graph-plotting-in-python-set-1/ inspired
#YOU HAVE TO HAVE INSTALLED THE FOLLOWING LIBRARIES : matplotlib.pyplot, colorama
# [pip install matplotlib.pyplot]
# [pip install colorama]
# normally alrady installed: math, time

import matplotlib.pyplot as plt
import math
from colorama import Fore
from colorama import Style
import time

### ANALYTE in ml ###
a = input(" Enter the constant volume of the analyte in ml: ")
print(f"{Fore.GREEN}The volume of the analyte is{Style.RESET_ALL}", a, f"{Fore.GREEN}ml{Style.RESET_ALL}")              
mol_a = input("Enter the molarity of the analyte: ")
print(f"{Fore.GREEN}The molarity of the analyte is{Style.RESET_ALL}", mol_a)
### TITRANT in ml ###
mol_t = input("Enter the molarity of the titrant: ")
t  = float(a)*float(mol_a) / float(mol_t)-0.00001
print(f"{Fore.GREEN}The molarity of the titrant is{Style.RESET_ALL}", mol_t)
# wait 1/2 a sec
print("calculating...")
time.sleep(.5)
print(f"{Fore.GREEN}The volume of the titrant under the equivalence point is{Style.RESET_ALL}",t, f"{Fore.GREEN}ml{Style.RESET_ALL}")
#Volumeof titrant in ml:
t0= 0.00001
t1= math.trunc(float(t)*10000)/10000
t2= float(t)*0.9999
t3= float(t)*0.999
t4= float(t)*0.99
t5= float(t)*0.97
t6= float(t)*0.9
t7= float(t)*0.7
t8= float(t)*0.5
t9= float(t)*0.3
t10= float(t)*0.1
t11= (float(t1))*1.0001
t12= (float(t1))*1.001
t13= (float(t1))*1.01
t14= (float(t1))*1.03
t15= (float(t1))*1.1
t16= (float(t1))*1.3
t17= (float(t1))*1.5
t18= (float(t1))*1.7
t19= (float(t1))*1.9
t20= (float(t1))*2
# mol of Analyte and Titrant:
anzahl_mol_von_starker_analyt = float(a)* 10**(-3)*float(mol_a)
#M                            =             z.B. 0.005      mol
anzahl_mol_von_starker_titrant = float(t)*10**(-3)*float(mol_t)
#pH at this ratio:
berechne_pH_t0       =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t0)*10**(-3)*float(mol_t)))       /       (float(t0)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_t1       =     ( 7)
berechne_pH_unter_t2 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t2)*10**(-3)*float(mol_t)))       /       (float(t2)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t3 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t3)*10**(-3)*float(mol_t)))       /       (float(t3)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t4 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t4)*10**(-3)*float(mol_t)))       /       (float(t4)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t5 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t5)*10**(-3)*float(mol_t)))       /       (float(t5)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t6 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t6)*10**(-3)*float(mol_t)))       /       (float(t6)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t7 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t7)*10**(-3)*float(mol_t)))       /       (float(t7)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t8 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t8)*10**(-3)*float(mol_t)))       /       (float(t8)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t9 =     (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t9)*10**(-3)*float(mol_t)))       /       (float(t9)*10**(-3)  + float(a)*10**(-3)   )       )
berechne_pH_unter_t10 =    (-1)*math.log10(        (float(anzahl_mol_von_starker_analyt) -    float(float(t10)*10**(-3)*float(mol_t)))      /       (float(t10)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t11 = 14-(-1)*math.log10(        (float(float(t11)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t11)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t12 = 14-(-1)*math.log10(        (float(float(t12)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t12)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t13 = 14-(-1)*math.log10(        (float(float(t13)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t13)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t14 = 14-(-1)*math.log10(        (float(float(t14)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t14)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t15 = 14-(-1)*math.log10(        (float(float(t15)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t15)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t16 = 14-(-1)*math.log10(        (float(float(t16)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t16)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t17 = 14-(-1)*math.log10(        (float(float(t17)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t17)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t18 = 14-(-1)*math.log10(        (float(float(t18)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t18)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t19 = 14-(-1)*math.log10(        (float(float(t19)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t19)*10**(-3) + float(a)*10**(-3)   )       )
berechne_pH_ueber_t20 = 14-(-1)*math.log10(        (float(float(t20)*10**(-3)*float(mol_t) - (float(anzahl_mol_von_starker_analyt))))       /       (float(t20)*10**(-3) + float(a)*10**(-3)   )       )
# X-axis: Volume of Titrant in ml:
x = [float(t0),float(t1),float(t2),float(t3),float(t4),float(t5),float(t6),float(t7),float(t8),float(t9),float(t10),float(t11),float(t12),float(t13),float(t14),float(t15),float(t16),float(t17),float(t18),float(t19),float(t20)]
# Y-axis: pH-value at this ratio:
y = [berechne_pH_t0,berechne_pH_t1,berechne_pH_unter_t2,berechne_pH_unter_t3,berechne_pH_unter_t4,berechne_pH_unter_t5,berechne_pH_unter_t6,berechne_pH_unter_t7,berechne_pH_unter_t8,berechne_pH_unter_t9,berechne_pH_unter_t10,berechne_pH_ueber_t11,berechne_pH_ueber_t12,berechne_pH_ueber_t13,berechne_pH_ueber_t14,berechne_pH_ueber_t15,berechne_pH_ueber_t16,berechne_pH_ueber_t17,berechne_pH_ueber_t18,berechne_pH_ueber_t19,berechne_pH_ueber_t20 ]
# Hiermit werden die Punkte in das Koordinatensystem geplottet:
plt.scatter(x, y, label= "values", color= "green", 
            marker= "o", s=30) 
# X-axis; title:
plt.xlabel('Volume of titrant in ml')
# Y-axis; title:
plt.ylabel('pH-value at this ratio')
# Title of graph:
plt.title('Titration curve, strong acid, strong base')
# show legend
plt.legend() 
# show plot
plt.show()

Solution

In order to fit the parameters of a parametric function to data, you can use scipy.optimize.curve_fit:

def my_function(x, A, B, C, D, E):
    return A/(1 + B**(x - C)) + D + E*x

parameters, covariance = curve_fit(f = my_function, xdata = x, ydata = y)

If necessary, you can improve the results by passing some bounds for the parameters:

parameters, covariance = curve_fit(f = my_function, xdata = x, ydata = y, bounds = ([0, 0, 0.9*t, -10, -10], [14, 1, 1.1*t, 10, 10]))

Complete Code

# IMPORT
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit


# INPUT VALUES
a = 100
mol_a = 10
mol_t = 5


# DATA GENERATION
N = 40
t  = a*mol_a/mol_t

x1 = np.linspace(0, t - 1e-6, N//2)
x2 = np.linspace(t + 1e-6, 2*t, N//2)
anzahl_mol_von_starker_analyt = a*10**(-3)*mol_a
anzahl_mol_von_starker_titrant = t*10**(-3)*mol_t

y1 = (-1)*np.log10((anzahl_mol_von_starker_analyt - x1*10**(-3)*mol_t)/(x1*10**(-3) + a*10**(-3)))
y2 = 14 - (-1)*np.log10((x2*10**(-3)*mol_t - anzahl_mol_von_starker_analyt)/(x2*10**(-3) + a*10**(-3)))

x = np.concatenate((x1, x2))
y = np.concatenate((y1, y2))


# PARAMETERS FITTING
def my_function(x, A, B, C, D, E):
    return A/(1 + B**(x - C)) + D + E*x

parameters, covariance = curve_fit(f = my_function, xdata = x, ydata = y, bounds = ([0, 0, 0.9*t, -10, -10], [14, 1, 1.1*t, 10, 10]))

for parameter, name in zip(parameters, ['A', 'B', 'C', 'D', 'E']):
    print(f'{name} = {parameter:14.10f}')

x_fitted = np.linspace(x[0], x[-1], 1000)
y_fitted = my_function(x_fitted, *parameters)


# PLOT
plt.style.use('seaborn-darkgrid')
fig, ax = plt.subplots()

ax.plot(x, y, label = 'data', marker = 'o', linestyle = '')
ax.plot(x_fitted, y_fitted, label = 'fitted curve')

ax.set_xlabel('Volume of titrant in ml')
ax.set_ylabel('pH-value at this ratio')
ax.set_title('Titration curve, strong acid, strong base')
ax.legend(frameon = True)

plt.show()

Result

A =  13.2787117583
B =   0.7311691018
C = 199.8189983959
D =  -0.8940392051
E =   0.0054250609

PLOT

Answered By - Zephyr