From hi, 3 Weeks ago, written in Python.
Embed
Hits: 108
1. import numpy as np
2. import matplotlib.pyplot as plt
3.
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_1 = np.zeros(n + 1)
15. a_1 = np.zeros(n)
16. b_1 = np.zeros(n)
17. c_1 = np.zeros(n)
18. d_1 = np.zeros(n)
19. M_2 = np.zeros(n + 1)
20. a_2 = np.zeros(n)
21. b_2 = np.zeros(n)
22. c_2 = np.zeros(n)
23. d_2 = np.zeros(n)
24.
25. z_matrix = z.reshape(-1, 1)
26.
27. matrix_1 = np.array([
28.     [4, 1, 0, 0, 0],
29.     [1, 4, 1, 0, 0],
30.     [0, 1, 4, 1, 0],
31.     [0, 0, 1, 4, 1],
32.     [0, 0, 0, 1, 4]
33. ])
34.
35. matrix_2 = np.array([
36.     [5, 1, 0, 0, 0],
37.     [1, 4, 1, 0, 0],
38.     [0, 1, 4, 1, 0],
39.     [0, 0, 1, 4, 1],
40.     [0, 0, 0, 1, 5]
41. ])
42.
43. inverse_matrix_1 = np.linalg.inv(matrix_1)
44. inverse_matrix_2 = np.linalg.inv(matrix_2)
45. M_matrix_1 = inverse_matrix_1 @ z_matrix
46. M_matrix_2 = inverse_matrix_2 @ z_matrix
47.
48. for i in range(n - 1):
49.     M_1[i + 1] = M_matrix_1[i].item()
50.     M_2[i + 1] = M_matrix_2[i].item()
51.
52. M_1[0] = 0
53. M_1[n] = 0
54. M_2[0] = M_2[1]
55. M_2[n] = M_2[n - 1]
56.
57. for i in range(n):
58.     a_1[i] = (M_1[i + 1] - M_1[i]) / (6 * h)
59.     b_1[i] = M_1[i] / 2
60.     c_1[i] = (y[i + 1] - y[i]) / h - h * (M_1[i + 1]
61.                                           + 2 * M_1[i]) / 6
62.     d_1[i] = y[i]
63.
64.     a_2[i] = (M_2[i + 1] - M_2[i]) / (6 * h)
65.     b_2[i] = M_2[i] / 2
66.     c_2[i] = (y[i + 1] - y[i]) / h - h * (M_2[i + 1]
67.                                           + 2 * M_2[i]) / 6
68.     d_2[i] = y[i]
69.
70. year_estimates = [1990, 2010]
71. spline_estimates_1 = np.zeros(2)
72. spline_estimates_2 = np.zeros(2)
73.
74. for j in range(2):
75.     year = year_estimates[j]
76.     i = min(n - 1, int((year - 1850) / h))
77.     dx = year - x[i]
78.     spline_estimates_1[j] = (a_1[i] * dx ** 3 + b_1[i] * dx ** 2
79.                              + c_1[i] * dx + d_1[i])
80.     spline_estimates_2[j] = (a_2[i] * dx ** 3 + b_2[i] * dx ** 2
81.                              + c_2[i] * dx + d_2[i])
82.
83. xx_extended = np.linspace(x[0], 2099, 1000)
84. yy_1_extended = np.zeros(len(xx_extended))
85. yy_2_extended = np.zeros(len(xx_extended))
86.
87. for k in range(len(xx_extended)):
88.     year = xx_extended[k]
89.     i = min(n - 1, int((year - 1850) / h))
90.     dx = year - x[i]
91.     yy_1_extended[k] = (a_1[i] * dx ** 3
92.                         + b_1[i] * dx ** 2
93.                         + c_1[i] * dx + d_1[i])
94.     yy_2_extended[k] = (a_2[i] * dx ** 3
95.                         + b_2[i] * dx ** 2
96.                         + c_2[i] * dx + d_2[i])
97.
98. x_1 = np.arange(1850, 2012)
99. y_1 = [
100.     285.2, 285.1, 285.0, 285.0, 284.9, 285.1, 285.4, 285.6,
101.     285.9, 286.1, 286.4, 286.6, 286.7, 286.8, 286.9, 287.1,
102.     287.2, 287.3, 287.4, 287.5, 287.7, 287.9, 288.0, 288.2,
103.     288.4, 288.6, 288.7, 288.9, 289.5, 290.1, 290.8, 291.4,
104.     292.0, 292.5, 292.9, 293.3, 293.8, 294.0, 294.1, 294.2,
105.     294.4, 294.6, 294.8, 294.7, 294.8, 294.8, 294.9, 294.9,
106.     294.9, 295.3, 295.7, 296.2, 296.6, 297.0, 297.5, 298.0,
107.     298.4, 298.8, 299.3, 299.7, 300.1, 300.6, 301.0, 301.3,
108.     301.4, 301.6, 302.0, 302.4, 302.8, 303.0, 303.4, 303.7,
109.     304.1, 304.5, 304.9, 305.3, 305.8, 306.2, 306.6, 307.2,
110.     307.5, 308.0, 308.3, 308.9, 309.3, 309.7, 310.1, 310.6,
111.     311.0, 311.2, 311.3, 311.0, 310.7, 310.5, 310.2, 310.3,
112.     310.3, 310.4, 310.5, 310.9, 311.3, 311.8, 312.2, 312.6,
113.     313.2, 313.7, 314.3, 314.8, 315.34, 316.18, 317.07,
114.     317.73, 318.43, 319.08, 319.65, 320.23, 321.59, 322.31,
115.     323.04, 324.23, 325.54, 326.42, 327.45, 329.43, 330.21,
116.     331.36, 331.92, 333.73, 335.42, 337.10, 338.99, 340.36,
117.     341.57, 342.53, 344.24, 345.72, 347.15, 348.93, 351.47,
118.     353.15, 354.29, 355.68, 356.42, 357.13, 358.61, 360.67,
119.     362.58, 363.48, 366.27, 368.38, 369.64, 371.15, 373.15,
120.     375.64, 377.44, 379.46, 381.59, 383.37, 385.46, 386.95,
121.     389.21, 391.15
122. ]
123.
124. x_2 = np.arange(2012, 2100)
125. y_2 = [
126.     389.9, 391.5, 393.1, 394.6, 396.2, 397.8, 399.4, 400.9,
127.     402.5, 404.0, 405.5, 407.0, 408.5, 410.0, 411.5, 413.0,
128.     414.5, 416.0, 417.5, 418.9, 420.4, 421.8, 423.2, 424.6,
129.     426.1, 427.5, 428.9, 430.2, 431.6, 433.0, 434.4, 435.7,
130.     437.1, 438.4, 439.8, 441.1, 442.4, 443.7, 445.0, 446.2,
131.     447.4, 448.5, 449.6, 450.8, 451.8, 452.9, 453.9, 454.9,
132.     455.9, 456.9, 457.8, 458.7, 459.6, 460.4, 461.3, 462.1,
133.     462.9, 463.6, 464.4, 465.1, 465.7, 466.4, 467.0, 467.6,
134.     468.2, 468.8, 469.3, 469.9, 470.4, 470.8, 471.2, 471.7,
135.     472.1, 472.4, 472.8, 473.1, 473.4, 473.6, 473.9, 474.1,
136.     474.3, 474.5, 474.6, 474.8, 474.9, 474.9, 475.0, 475.0
137. ]
138.
139. x_3 = np.arange(2012, 2100)
140. y_3 = {
141.     392.0, 394.0, 396.0, 398.0, 400.1, 402.3, 404.4, 406.5,
142.     408.7, 410.9, 413.1, 415.4, 417.7, 420.0, 422.4, 424.7,
143.     427.1, 429.5, 432.0, 434.4, 436.9, 439.5, 442.0, 444.5,
144.     447.1, 449.8, 452.5, 455.1, 457.8, 460.5, 463.2, 466.0,
145.     468.9, 471.7, 474.5, 477.4, 480.3, 483.3, 486.2, 489.2,
146.     492.1, 494.9, 497.6, 500.4, 503.0, 505.6, 508.1, 510.6,
147.     513.0, 515.4, 517.7, 519.8, 522.0, 524.1, 526.2, 528.1,
148.     530.1, 532.0, 533.8, 535.5, 537.2, 538.8, 540.3, 541.8,
149.     543.3, 544.7, 546.0, 547.3, 548.5, 549.7, 550.8, 551.7,
150.     552.7, 553.6, 554.5, 555.2, 556.0, 556.7, 557.2, 557.8,
151.     558.3, 558.7, 559.0, 559.3, 559.6, 559.8, 560.0, 560.0
152. }
153.
154. plt.figure(figsize=(12, 6))
155. plt.plot(x, y, 'ro', label='Data points')
156. plt.plot(xx_extended, yy_1_extended, 'r-',
157.          label='Natural Spline')
158. plt.plot(xx_extended, yy_2_extended, 'g-',
159.          label='Parabolic Runout Spline')
160. plt.plot(x_1, y_1, 'b-', label='Real Data')
161. plt.plot(x_2, y_2, 'm-', label='Future Scenario 1')
162. plt.plot(x_3, y_3, 'c-', label='Future Scenario 2')
163.
164. plt.xlabel('Year')
165. plt.ylabel('CO2 Concentration (ppm)')
166. plt.legend(loc='best')
167. plt.grid(True)
168. plt.title('Global Mean CO2 Mixing Ratio over Time')
169. plt.show()
170.