fitting theory with experimental data

I want to do fitting on experimental data and my theory I want to try to find the parameters that fit the same as the theory and experiment do I have to make an algorithm? what algorithm do I use to get a good fitting and mse for each parameter which is very small

I avoid using modules to do my fitting as much as possible trying to understand how the code works.

I use this mse

    N = len(y)
error = (np.sum(y**2) - 2 * np.sum(y * y_t_fit) + np.sum(y_t_fit**2)) / N

this my code

import numpy as np
import matplotlib.pyplot as plt


#data
time = np.array([
    0.000,
    0.034,
    0.068,
    0.101,
    0.135,
    0.169,
    0.203,
    0.236,
    0.270,
    0.304,
    0.338,
    0.372,
    0.405,
    0.439,
    0.473,
    0.507,
    0.541,
    0.574,
    0.608,
    0.642,
    0.676,
    0.709,
    0.743,
    0.777,
    0.811,
    0.845,
    0.878,
    0.912,
    0.946,
    0.980,
    1.013,
    1.047,
    1.081,
    1.115,
    1.149,
    1.182,
    1.216,
    1.250,
    1.284,
    1.317,
    1.351,
    1.385,
    1.419,
    1.453,
    1.486,
    1.520,
    1.554,
    1.588,
    1.622,
    1.655,
    1.689,
    1.723,
    1.757,
    1.790,
    1.824,
    1.858,
    1.892,
    1.926,
    1.959,
    1.993,
    2.027,
    2.061,
    2.094,
    2.128,
    2.162,
    2.196,
    2.230,
    2.263,
    2.297,
    2.331,
    2.365,
    2.398,
    2.432,
    2.466,
    2.500,
    2.534,
    2.567,
    2.601,
    2.635,
    2.669,
    2.703,
    2.736,
    2.770,
    2.804,
    2.838,
    2.871,
    2.905,
    2.939,
    2.973,
    3.007,
    3.040,
    3.074,
    3.108,
    3.142,
    3.145
])

y = np.array([
    0.827,
    0.828,
    0.827,
    0.817,
    0.798,
    0.769,
    0.731,
    0.692,
    0.651,
    0.607,
    0.561,
    0.510,
    0.469,
    0.435,
    0.397,
    0.368,
    0.344,
    0.325,
    0.316,
    0.314,
    0.324,
    0.333,
    0.346,
    0.372,
    0.399,
    0.429,
    0.457,
    0.485,
    0.519,
    0.545,
    0.566,
    0.585,
    0.602,
    0.618,
    0.628,
    0.650,
    0.636,
    0.619,
    0.608,
    0.592,
    0.576,
    0.559,
    0.538,
    0.513,
    0.493,
    0.478,
    0.460,
    0.444,
    0.429,
    0.422,
    0.415,
    0.413,
    0.413,
    0.417,
    0.430,
    0.441,
    0.450,
    0.465,
    0.478,
    0.493,
    0.506,
    0.518,
    0.529,
    0.536,
    0.543,
    0.550,
    0.548,
    0.546,
    0.544,
    0.539,
    0.530,
    0.520,
    0.510,
    0.504,
    0.496,
    0.486,
    0.476,
    0.471,
    0.469,
    0.463,
    0.457,
    0.456,
    0.458,
    0.459,
    0.464,
    0.471,
    0.479,
    0.484,
    0.489,
    0.497,
    0.505,
    0.507,
    0.511,
    0.517,
    0.516
])

y -= 0.5

# find the max
def max_p (y):
    result = []
    for i in range(1,len(y)-1):
        if y[i-1] < y[i] > y[i+1]:
            result.append(i)
    return result
max = max_p(y)

#T1
y0 = y[max[0]]
yT = y[max[1]]
T = time[max[1]] - time[max[0]] 

#T2
y1 = y[max[1]]
yT1 = y[max[2]]
T = time[max[2]] - time[max[1]]

max = max_p(y)

peak_values = [y[i] for i in max]



#Parameter
A = 0.32
b = 0.5
k = 4
m = 1

w = np.sqrt(k/m - (b/2*m)**2)
T = 2*np.pi/w

y_t = np.exp((-b/(2*m))*t) * (A*np.cos(w*t))


plt.plot(time,y,label='experiment',linestyle='-')
plt.scatter(time[max], peak_values, color='red', marker='o', s=100, label="Max")
plt.plot(t,y_t,"r--")
plt.legend()
plt.grid()
plt.show()

You are missing scipy in your Python process. Lets import libraries and your data.

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

t = np.array([ 0.000, 0.034, 0.068, 0.101, 0.135, 0.169, 0.203, 0.236, 0.270, 0.304, 0.338, 0.372, 0.405, 0.439, 0.473, 0.507, 0.541, 0.574, 0.608, 0.642, 0.676, 0.709, 0.743, 0.777, 0.811, 0.845, 0.878, 0.912, 0.946, 0.980, 1.013, 1.047, 1.081, 1.115, 1.149, 1.182, 1.216, 1.250, 1.284, 1.317, 1.351, 1.385, 1.419, 1.453, 1.486, 1.520, 1.554, 1.588, 1.622, 1.655, 1.689, 1.723, 1.757, 1.790, 1.824, 1.858, 1.892, 1.926, 1.959, 1.993, 2.027, 2.061, 2.094, 2.128, 2.162, 2.196, 2.230, 2.263, 2.297, 2.331, 2.365, 2.398, 2.432, 2.466, 2.500, 2.534, 2.567, 2.601, 2.635, 2.669, 2.703, 2.736, 2.770, 2.804, 2.838, 2.871, 2.905, 2.939, 2.973, 3.007, 3.040, 3.074, 3.108, 3.142, 3.145 ])
y = np.array([ 0.827, 0.828, 0.827, 0.817, 0.798, 0.769, 0.731, 0.692, 0.651, 0.607, 0.561, 0.510, 0.469, 0.435, 0.397, 0.368, 0.344, 0.325, 0.316, 0.314, 0.324, 0.333, 0.346, 0.372, 0.399, 0.429, 0.457, 0.485, 0.519, 0.545, 0.566, 0.585, 0.602, 0.618, 0.628, 0.650, 0.636, 0.619, 0.608, 0.592, 0.576, 0.559, 0.538, 0.513, 0.493, 0.478, 0.460, 0.444, 0.429, 0.422, 0.415, 0.413, 0.413, 0.417, 0.430, 0.441, 0.450, 0.465, 0.478, 0.493, 0.506, 0.518, 0.529, 0.536, 0.543, 0.550, 0.548, 0.546, 0.544, 0.539, 0.530, 0.520, 0.510, 0.504, 0.496, 0.486, 0.476, 0.471, 0.469, 0.463, 0.457, 0.456, 0.458, 0.459, 0.464, 0.471, 0.479, 0.484, 0.489, 0.497, 0.505, 0.507, 0.511, 0.517, 0.516 ])

Then we define the following model:

def model(t, a, rho, p, omega, phi):
    return a + rho * np.exp(- p * t) * np.sin(2 * np.pi * omega * t + phi)

And we perform the optimization:

popt, pcov = optimize.curve_fit(model, t, y)
# array([0.5006398 , 0.34257547, 0.85005221, 0.91936647, 0.96196433])

Which renders as follow:

tlin = np.linspace(t.min(), t.max(), 200)
ylin = model(tlin, *popt)

As you can see, the general shape of the dataset is captured but we cannot say that it is a perfect match, probably because your dataset is not exactly compliant with the model.