1. import numpy as np
2. import matplotlib.pyplot as plt
3. # a
4. x = np.array([1850, 1875, 1900, 1925, 1950, 1975, 2000])
5. y = np.array([285.2, 288.6, 295.7, 305.3, 311.3, 331.36,
6.               369.64])
7. h = 25
8. n = len(x) - 1
9.  
10. z = np.zeros(n-1)
11. for i in range(n-1):
12.     z[i] = (y[i] - 2 * y[i+1] + y[i+2]) * 6 / (h**2)
13.  
14. M = np.zeros(n+1)
15. a = np.zeros(n)
16. b = np.zeros(n)
17. c = np.zeros(n)
18. d = np.zeros(n)
19.  
20. z_matrix = z
21.  
22. matrix = (np.diag(4*np.ones(n-1)) + np.diag(np.ones(n-2), 1)
23.           + np.diag(np.ones(n-2), -1))
24.  
25. inverse_matrix = np.linalg.inv(matrix)
26. M_matrix = np.dot(inverse_matrix, z_matrix)
27. for i in range(n-1):
28.     M[i+1] = M_matrix[i]
29. M[0] = 0
30. M[n] = 0
31.  
32. for i in range(n):
33.     a[i] = (M[i+1] - M[i]) / (6 * h)
34.     b[i] = M[i] / 2
35.     c[i] = (y[i+1] - y[i]) / h - h * (M[i+1] + 2 * M[i]) / 6
36.     d[i] = y[i]
37.  
38. for i in range(n):
39.     print(f"S_{i+1}(x) = {a[i]:.4e}(x - {x[i]})^3 "
40.           f"+ {b[i]:.4e}(x - {x[i]})^2 + {c[i]:.4e}(x - {x[i]}) "
41.           f"+ {d[i]:.4e}, for {x[i]} <= x < {x[i+1]}")
42.  
43. year_estimates = [1990, 2010]
44. spline_estimates = np.zeros(2)
45. for j in range(2):
46.     year = year_estimates[j]
47.     i = int((year - 1850) / h)
48.     if i >= n:
49.         i = n - 1
50.     dx = year - x[i]
51.     spline_estimates[j] = (a[i]*dx**3 + b[i]*dx**2
52.                            + c[i]*dx + d[i])
53.  
54. print(f"\nCO2 concentration estimate for 1990="
55.       f" {spline_estimates[0]:.2f} ppm")
56. print(f"CO2 concentration estimate for 2010= "
57.       f"{spline_estimates[1]:.2f} ppm")
58.  
59. print("a coefficients=", a)
60. print("b coefficients=", b)
61. print("c coefficients=", c)
62. print("d coefficients=", d)
63. print("M values=", M)
64.  
65. xx = np.linspace(x[0], x[-1], 1000)
66. yy = np.zeros(len(xx))
67. for k in range(len(xx)):
68.     year = xx[k]
69.     i = int((year - 1850) / h)
70.     if i >= n:
71.         i = n - 1
72.     dx = year - x[i]
73.     yy[k] = a[i]*dx**3 + b[i]*dx**2 + c[i]*dx + d[i]
74.  
75. plt.figure(figsize=(12, 6))
76. plt.plot(x, y, 'o', label='Original data')
77. plt.plot(xx, yy, label='Cubic spline')
78. plt.xlabel('Year')
79. plt.ylabel('CO2 Concentration (ppm)')
80. plt.title('Natural Cubic Spline Interpolation of CO2 Concentration')
81. plt.legend()
82. plt.grid(True)
83. plt.show()
84.