How to do fitting on experimental graph and theoretical graph

experimental and theoretical graphs

I am currently fitting the experimental results and will later take the MSE to increase the accuracy, the accuracy between the data and theory.

How do you make sure the fittings fit perfectly to each other? Can that be done? And how do you calculate the MSE?

so this is my code

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


A = 0.12
b = 0.005     
m = 1.0      
k = 7.0      



gamma = b / (2 * m)

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.175, 3.209, 3.243
])

y = np.array([
    0.309, 0.320, 0.326, 0.325, 0.321, 0.312, 0.299, 0.282, 0.266, 0.249,
    0.230, 0.209, 0.187, 0.169, 0.154, 0.138, 0.126, 0.115, 0.107, 0.103,
    0.102, 0.107, 0.111, 0.116, 0.128, 0.140, 0.153, 0.165, 0.177, 0.192,
    0.203, 0.212, 0.221, 0.228, 0.235, 0.239, 0.239, 0.240, 0.236, 0.231,
    0.224, 0.217, 0.210, 0.200, 0.189, 0.181, 0.174, 0.166, 0.159, 0.152,
    0.149, 0.146, 0.145, 0.145, 0.147, 0.153, 0.157, 0.162, 0.167, 0.173,
    0.180, 0.186, 0.191, 0.196, 0.199, 0.202, 0.205, 0.205, 0.204, 0.203,
    0.201, 0.197, 0.193, 0.189, 0.186, 0.182, 0.178, 0.174, 0.171, 0.171,
    0.168, 0.166, 0.165, 0.166, 0.166, 0.169, 0.172, 0.175, 0.177, 0.179,
    0.182, 0.186, 0.186, 0.188, 0.190, 0.192, 0.192
])

y -= 0.2  # Koreksi offset

peaks, _ = find_peaks(y)
periode = np.diff(time[peaks])



print("Indeks puncak:", peaks)
print("Waktu puncak:", time[peaks])
print("Periode antara puncak:", periode)
print("Periode rata-rata:", np.mean(periode))

# Ambil puncak pertama dan terakhir
y0 = y[peaks[0]]
yT = y[peaks[-1]]
T = time[peaks[-1]] - time[peaks[0]]  # Selisih waktu antara puncak pertama dan terakhir


# Mencari puncak
peaks, _ = find_peaks(y)
periode = np.diff(time[peaks])

# Hitung omega_d berdasarkan data
omega_d = 2 * np.pi / np.mean(periode)

# Fungsi fitting dengan gamma sebagai parameter
def damped_cosine(t, A_fit, gamma_fit, phi_fit):
    return A_fit * np.exp(-gamma_fit * t) * np.cos(omega_d * t + phi_fit)

# Fitting ke data eksperimen
popt, _ = curve_fit(damped_cosine, time, y, p0=[0.12, gamma, 0])

# Ambil hasil fitting
A_fit, gamma_fit, phi_fit = popt

# Hitung ulang data teori dengan hasil fitting
y_t_fit = damped_cosine(time, A_fit, gamma_fit, phi_fit)

# Plot hasil
plt.plot(time, y, label="Experiment")
plt.plot(time[peaks], y[peaks], "ro", label="Max")
plt.plot(time, y_t_fit, label="Fitting", linestyle="dashed")
plt.xlabel("time (s)")
plt.ylabel(" y(t)")
plt.title("Grafik Getaran Teredam dengan Fitting")
plt.legend()
plt.grid()
plt.show()

# Cetak hasil fitting
print(f"Amplitudo fit: {A_fit}")
print(f"Gamma fit: {gamma_fit}")
print(f"Fase fit: {phi_fit}")


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

# Cetak error
print(f"Error rata-rata kuadratik: {error}")

Solution

Your experimental data has a negative offset in it (not centred on y=0). You should allow for that in your fitting function. Although the peaks are a good first indicator of period, the damped frequency will differ from the undamped one, so I would allow that to vary also.

Your error estimation looks OK (if a little long).

Changed lines (with full code below):

def damped_cosine(t, z_fit, A_fit, omega_fit, gamma_fit, phi_fit):
    return z_fit + A_fit * np.exp(-gamma_fit * t) * np.cos(omega_fit * t + phi_fit)

popt, _ = curve_fit(damped_cosine, time, y, p0=[0.0, 0.12, omega_d, gamma, 0])
z_fit, A_fit, omega_fit, gamma_fit, phi_fit = popt

y_t_fit = damped_cosine(time, z_fit, A_fit, omega_fit, gamma_fit, phi_fit)

Full code:

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


A = 0.12
b = 0.005     
m = 1.0      
k = 7.0      



gamma = b / (2 * m)

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.175, 3.209, 3.243
])

y = np.array([
    0.309, 0.320, 0.326, 0.325, 0.321, 0.312, 0.299, 0.282, 0.266, 0.249,
    0.230, 0.209, 0.187, 0.169, 0.154, 0.138, 0.126, 0.115, 0.107, 0.103,
    0.102, 0.107, 0.111, 0.116, 0.128, 0.140, 0.153, 0.165, 0.177, 0.192,
    0.203, 0.212, 0.221, 0.228, 0.235, 0.239, 0.239, 0.240, 0.236, 0.231,
    0.224, 0.217, 0.210, 0.200, 0.189, 0.181, 0.174, 0.166, 0.159, 0.152,
    0.149, 0.146, 0.145, 0.145, 0.147, 0.153, 0.157, 0.162, 0.167, 0.173,
    0.180, 0.186, 0.191, 0.196, 0.199, 0.202, 0.205, 0.205, 0.204, 0.203,
    0.201, 0.197, 0.193, 0.189, 0.186, 0.182, 0.178, 0.174, 0.171, 0.171,
    0.168, 0.166, 0.165, 0.166, 0.166, 0.169, 0.172, 0.175, 0.177, 0.179,
    0.182, 0.186, 0.186, 0.188, 0.190, 0.192, 0.192
])

y -= 0.2  # Koreksi offset

peaks, _ = find_peaks(y)
periode = np.diff(time[peaks])



print("Indeks puncak:", peaks)
print("Waktu puncak:", time[peaks])
print("Periode antara puncak:", periode)
print("Periode rata-rata:", np.mean(periode))

# Ambil puncak pertama dan terakhir
y0 = y[peaks[0]]
yT = y[peaks[-1]]
T = time[peaks[-1]] - time[peaks[0]]  # Selisih waktu antara puncak pertama dan terakhir


# Mencari puncak
peaks, _ = find_peaks(y)
periode = np.diff(time[peaks])

# Hitung omega_d berdasarkan data
omega_d = 2 * np.pi / np.mean(periode)

# Fungsi fitting dengan gamma sebagai parameter
def damped_cosine(t, z_fit, A_fit, omega_fit, gamma_fit, phi_fit):
    return z_fit + A_fit * np.exp(-gamma_fit * t) * np.cos(omega_fit * t + phi_fit)

# Fitting ke data eksperimen
popt, _ = curve_fit(damped_cosine, time, y, p0=[0.0, 0.12, omega_d, gamma, 0])

# Ambil hasil fitting
z_fit, A_fit, omega_fit, gamma_fit, phi_fit = popt

# Hitung ulang data teori dengan hasil fitting
y_t_fit = damped_cosine(time, z_fit, A_fit, omega_fit, gamma_fit, phi_fit)

# Plot hasil
plt.plot(time, y, label="Experiment")
plt.plot(time[peaks], y[peaks], "ro", label="Max")
plt.plot(time, y_t_fit, label="Fitting", linestyle="dashed")
plt.xlabel("time (s)")
plt.ylabel(" y(t)")
plt.title("Grafik Getaran Teredam dengan Fitting")
plt.legend()
plt.grid()
plt.show()

# Cetak hasil fitting
print(f"Amplitudo fit: {A_fit}")
print(f"Gamma fit: {gamma_fit}")
print(f"Fase fit: {phi_fit}")


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

# Cetak error
print(f"Error rata-rata kuadratik: {error}")