Midterm Exam #1 (20232) Hands On Session Solutions

FIZ228 - Numerical Analysis
Dr. Emre S. Tasci, Hacettepe University

19/04/2024

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

1

Import the given FIZ228_20232_MT1_data.csv file as a pandas dataframe.

This data is known to contain signals taken from two joint sinusoidal signal generators (where superposition applies, i.e., $y = y_1 + y_2$).

Rename the first column as "t" and the second as "y"

data = pd.read_csv("../data/FIZ228_20232_MT1_data.csv",header=None)
data.columns = ["t","y"]
data

	t	y
0	0.000000	-0.631523
1	0.033445	-0.447049
2	0.066890	-0.331072
3	0.100334	-0.114149
4	0.133779	-0.168711
...	...	...
295	9.866221	-1.248320
296	9.899666	-0.943883
297	9.933110	-1.443523
298	9.966555	-1.640351
299	10.000000	-1.792804

300 rows × 2 columns

2

Plot its graph as using red dots (no lines).

plt.plot(data["t"],data["y"],"r.")
plt.show()

3

As mentioned above, this data consists of the superposition of two sinusoidal signals, such that:

$$y = A_1\sin(\omega_1 t + \varphi_1)+A_2\sin(\omega_2 t + \varphi_2)$$

Determine the parameters (amplitude, angular frequency and phase) of these two sinusoidal signals that best match to the data given.

(Hint: Depending on your method, you can relax the default tolerance, in other words, don't get upset if a very high precision match can't be obtained)

def f(t,params):
    (A1,A2,w1,w2,ph1,ph2) = params
    return A1*np.sin(w1*t+np.deg2rad(ph1)) + \
           A2*np.sin(w2*t+np.deg2rad(ph2))

def err_f(params):
    return ((data["y"] - f(data["t"],params))**2).sum()

res = opt.minimize(err_f,[3,2,2,1,0,-1])
res

      fun: 9.212506361972885
 hess_inv: array([[ 3.57003922e-03,  2.25792574e-04,  2.82854645e-05,
         1.33087083e-04, -1.45580153e-02, -3.44404868e-02],
       [ 2.25792574e-04,  4.49656576e-03, -1.38716170e-04,
        -4.40773031e-05,  7.99209961e-03,  2.34885404e-02],
       [ 2.82854645e-05, -1.38716170e-04,  6.20287290e-05,
         2.78434072e-05, -4.56936653e-03, -1.11882206e-02],
       [ 1.33087083e-04, -4.40773031e-05,  2.78434072e-05,
         8.61266226e-05, -9.93953936e-03, -2.47799129e-02],
       [-1.45580153e-02,  7.99209961e-03, -4.56936653e-03,
        -9.93953936e-03,  1.60636606e+00,  3.89858320e+00],
       [-3.44404868e-02,  2.34885404e-02, -1.11882206e-02,
        -2.47799129e-02,  3.89858320e+00,  9.74387856e+00]])
      jac: array([-4.76837158e-07, -3.09944153e-06,  2.24113464e-05,  1.39474869e-05,
        3.57627869e-07, -2.38418579e-07])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 371
      nit: 27
     njev: 53
   status: 2
  success: False
        x: array([  1.50023154,   2.01452057,   2.30302471,   0.99961517,
       -35.64788195,  -0.05551862])

4

Using these parameters, construct your estimation function and introduce a third column labeled "e" on your dataframe that holds the function values with respect to your t values.

data["e"] = f(data["t"],res.x)

data

	t	y	e
0	0.000000	-0.631523	-0.876290
1	0.033445	-0.447049	-0.712552
2	0.066890	-0.331072	-0.544274
3	0.100334	-0.114149	-0.372130
4	0.133779	-0.168711	-0.196818
...	...	...	...
295	9.866221	-1.248320	-1.014949
296	9.899666	-0.943883	-1.189764
297	9.933110	-1.443523	-1.361917
298	9.966555	-1.640351	-1.530671
299	10.000000	-1.792804	-1.695307

300 rows × 3 columns

5

Plot the given data against the model estimation (data points: red dots as above; model estimation: dashed line)

plt.plot(data["t"],data["e"],"b--")
plt.plot(data["t"],data["y"],"r.")
plt.show()

6

Calculate the following error estimations between the given data and your model:$\DeclareMathOperator\erf{erf}$

• $\erf_1 = \sqrt{\sum_{i}{(y_i - e_i)^2}}$

• $\erf_2 = \sum_{i}{|y_i - e_i|}$

  where $y$ indicates the given data and $e$ the model estimation.

  Also calculate the coefficient of determination ($r^2$).

erf1 = np.sqrt(((data["y"] - data["e"])**2).sum())
erf2 = np.abs((data["y"] - data["e"]).sum())
erf1,erf2

(3.0352110901834957, 3.886792829550611)

S_r = ((data["y"] - data["e"])**2).sum()
S_t = ((data["y"] - data["y"].mean())**2).sum()
r2 = (S_t - S_r) / S_t

r2

0.9890390943345185

7

The data presented was obtained by adding noise to the data calculated with the following parameters:

$$\begin{align*}A_1 &= 2.0& \omega_1 &= 1.0& \varphi_1 &=0.0\\ A_2 &= 1.5& \omega_2 &= 2.3& \varphi_2 &=-35.0^o \end{align*}$$

(You can't use these values to spot a good starting position for the above operations!)

Plot the data along with your fit and the model constructed with these parameters like you did in part 5.

Calculate the error estimations between the data and a model constructed with these parameters like you did in part 6.

params_org = [2.0,1.5,1.0,2.3,0.0,-35.0]
y_org = f(data["t"],params_org)

plt.plot(data["t"],y_org,"k-")
plt.plot(data["t"],data["e"],"b--")
plt.plot(data["t"],data["y"],"r.")
plt.legend(["original","fit","data"])
plt.show()

erf1o = np.sqrt(((data["y"] - y_org)**2).sum())
erf2o = np.abs((data["y"] - y_org).sum())
erf1o,erf2o

(3.0482001175866738, 3.222162372189898)

S_ro = ((data["y"] - y_org)**2).sum()
S_to = ((data["y"] - data["y"].mean())**2).sum()
r2o = (S_to - S_ro) / S_to

r2o

0.9889450803529385