1. import numpy as np
2. import matplotlib.pyplot as plt
3. # b
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.array([
23.     [5, 1, 0, 0, 0],
24.     [1, 4, 1, 0, 0],
25.     [0, 1, 4, 1, 0],
26.     [0, 0, 1, 4, 1],
27.     [0, 0, 0, 1, 5]
28. ])
29.  
30. inverse_matrix = np.linalg.inv(matrix)
31. M_matrix = np.dot(inverse_matrix, z_matrix)
32. for i in range(n-1):
33.     M[i+1] = M_matrix[i]
34. M[0] = M[1]
35. M[n] = M[n-1]
36.  
37. for i in range(n):
38.     a[i] = (M[i+1] - M[i]) / (6 * h)
39.     b[i] = M[i] / 2
40.     c[i] = (y[i+1] - y[i]) / h - h * (M[i+1] + 2 * M[i]) / 6
41.     d[i] = y[i]
42.  
43. for i in range(n):
44.     print(f"S_{i+1}(x) = {a[i]:.4e}(x - {x[i]})^3 "
45.           f"+ {b[i]:.4e}(x - {x[i]})^2 + {c[i]:.4e}(x - {x[i]}) "
46.           f"+ {d[i]:.4e}, for {x[i]} <= x < {x[i+1]}")
47.  
48. year_estimates = [1990, 2010]
49. spline_estimates = np.zeros(2)
50. for j in range(2):
51.     year = year_estimates[j]
52.     i = int((year - 1850) / h)
53.     if i >= n:
54.         i = n - 1
55.     dx = year - x[i]
56.     spline_estimates[j] = (a[i]*dx**3 + b[i]*dx**2
57.                            + c[i]*dx + d[i])
58.  
59. print(f"\nCO2 concentration estimate for 1990: "
60.       f"{spline_estimates[0]:.2f} ppm")
61. print(f"CO2 concentration estimate for 2010: "
62.       f"{spline_estimates[1]:.2f} ppm")
63.  
64. print("a coefficients=", a)
65. print("b coefficients=", b)
66. print("c coefficients=", c)
67. print("d coefficients=", d)
68. print("M values=", M)
69.  
70. xx = np.linspace(x[0], x[-1], 1000)
71. yy = np.zeros(len(xx))
72. for k in range(len(xx)):
73.     year = xx[k]
74.     i = int((year - 1850) / h)
75.     if i >= n:
76.         i = n - 1
77.     dx = year - x[i]
78.     yy[k] = a[i]*dx**3 + b[i]*dx**2 + c[i]*dx + d[i]
79.  
80. plt.figure(figsize=(12, 6))
81. plt.plot(x, y, 'o', label='Original data')
82. plt.plot(xx, yy, label='Parabolic Runout Cubic Spline')
83. plt.xlabel('Year')
84. plt.ylabel('CO2 Concentration (ppm)')
85. plt.title('Parabolic Runout Cubic Spline Interpolation of CO2 Concentration')
86. plt.legend()
87. plt.grid(True)
88. plt.show()
89.