import numpy as np import matplotlib.pyplot as plt x = np.array([1850, 1875, 1900, 1925, 1950, 1975, 2000]) y = np.array([285.2, 288.6, 295.7, 305.3, 311.3, 331.36, 369.64]) h = 25 n = len(x) - 1 z = np.zeros(n - 1) for i in range(n - 1): z[i] = (y[i] - 2 * y[i + 1] + y[i + 2]) * 6 / (h ** 2) M_1 = np.zeros(n + 1) a_1 = np.zeros(n) b_1 = np.zeros(n) c_1 = np.zeros(n) d_1 = np.zeros(n) M_2 = np.zeros(n + 1) a_2 = np.zeros(n) b_2 = np.zeros(n) c_2 = np.zeros(n) d_2 = np.zeros(n) z_matrix = z.reshape(-1, 1) matrix_1 = np.array([ [4, 1, 0, 0, 0], [1, 4, 1, 0, 0], [0, 1, 4, 1, 0], [0, 0, 1, 4, 1], [0, 0, 0, 1, 4] ]) matrix_2 = np.array([ [5, 1, 0, 0, 0], [1, 4, 1, 0, 0], [0, 1, 4, 1, 0], [0, 0, 1, 4, 1], [0, 0, 0, 1, 5] ]) inverse_matrix_1 = np.linalg.inv(matrix_1) inverse_matrix_2 = np.linalg.inv(matrix_2) M_matrix_1 = inverse_matrix_1 @ z_matrix M_matrix_2 = inverse_matrix_2 @ z_matrix for i in range(n - 1): M_1[i + 1] = M_matrix_1[i].item() M_2[i + 1] = M_matrix_2[i].item() M_1[0] = 0 M_1[n] = 0 M_2[0] = M_2[1] M_2[n] = M_2[n - 1] for i in range(n): a_1[i] = (M_1[i + 1] - M_1[i]) / (6 * h) b_1[i] = M_1[i] / 2 c_1[i] = (y[i + 1] - y[i]) / h - h * (M_1[i + 1] + 2 * M_1[i]) / 6 d_1[i] = y[i] a_2[i] = (M_2[i + 1] - M_2[i]) / (6 * h) b_2[i] = M_2[i] / 2 c_2[i] = (y[i + 1] - y[i]) / h - h * (M_2[i + 1] + 2 * M_2[i]) / 6 d_2[i] = y[i] year_estimates = [1990, 2010] spline_estimates_1 = np.zeros(2) spline_estimates_2 = np.zeros(2) for j in range(2): year = year_estimates[j] i = min(n - 1, int((year - 1850) / h)) dx = year - x[i] spline_estimates_1[j] = (a_1[i] * dx ** 3 + b_1[i] * dx ** 2 + c_1[i] * dx + d_1[i]) spline_estimates_2[j] = (a_2[i] * dx ** 3 + b_2[i] * dx ** 2 + c_2[i] * dx + d_2[i]) xx_extended = np.linspace(x[0], 2099, 1000) yy_1_extended = np.zeros(len(xx_extended)) yy_2_extended = np.zeros(len(xx_extended)) for k in range(len(xx_extended)): year = xx_extended[k] i = min(n - 1, int((year - 1850) / h)) dx = year - x[i] yy_1_extended[k] = (a_1[i] * dx ** 3 + b_1[i] * dx ** 2 + c_1[i] * dx + d_1[i]) yy_2_extended[k] = (a_2[i] * dx ** 3 + b_2[i] * dx ** 2 + c_2[i] * dx + d_2[i]) x_1 = np.arange(1850, 2012) y_1 = [ 285.2, 285.1, 285.0, 285.0, 284.9, 285.1, 285.4, 285.6, 285.9, 286.1, 286.4, 286.6, 286.7, 286.8, 286.9, 287.1, 287.2, 287.3, 287.4, 287.5, 287.7, 287.9, 288.0, 288.2, 288.4, 288.6, 288.7, 288.9, 289.5, 290.1, 290.8, 291.4, 292.0, 292.5, 292.9, 293.3, 293.8, 294.0, 294.1, 294.2, 294.4, 294.6, 294.8, 294.7, 294.8, 294.8, 294.9, 294.9, 294.9, 295.3, 295.7, 296.2, 296.6, 297.0, 297.5, 298.0, 298.4, 298.8, 299.3, 299.7, 300.1, 300.6, 301.0, 301.3, 301.4, 301.6, 302.0, 302.4, 302.8, 303.0, 303.4, 303.7, 304.1, 304.5, 304.9, 305.3, 305.8, 306.2, 306.6, 307.2, 307.5, 308.0, 308.3, 308.9, 309.3, 309.7, 310.1, 310.6, 311.0, 311.2, 311.3, 311.0, 310.7, 310.5, 310.2, 310.3, 310.3, 310.4, 310.5, 310.9, 311.3, 311.8, 312.2, 312.6, 313.2, 313.7, 314.3, 314.8, 315.34, 316.18, 317.07, 317.73, 318.43, 319.08, 319.65, 320.23, 321.59, 322.31, 323.04, 324.23, 325.54, 326.42, 327.45, 329.43, 330.21, 331.36, 331.92, 333.73, 335.42, 337.10, 338.99, 340.36, 341.57, 342.53, 344.24, 345.72, 347.15, 348.93, 351.47, 353.15, 354.29, 355.68, 356.42, 357.13, 358.61, 360.67, 362.58, 363.48, 366.27, 368.38, 369.64, 371.15, 373.15, 375.64, 377.44, 379.46, 381.59, 383.37, 385.46, 386.95, 389.21, 391.15 ] x_2 = np.arange(2012, 2100) y_2 = [ 389.9, 391.5, 393.1, 394.6, 396.2, 397.8, 399.4, 400.9, 402.5, 404.0, 405.5, 407.0, 408.5, 410.0, 411.5, 413.0, 414.5, 416.0, 417.5, 418.9, 420.4, 421.8, 423.2, 424.6, 426.1, 427.5, 428.9, 430.2, 431.6, 433.0, 434.4, 435.7, 437.1, 438.4, 439.8, 441.1, 442.4, 443.7, 445.0, 446.2, 447.4, 448.5, 449.6, 450.8, 451.8, 452.9, 453.9, 454.9, 455.9, 456.9, 457.8, 458.7, 459.6, 460.4, 461.3, 462.1, 462.9, 463.6, 464.4, 465.1, 465.7, 466.4, 467.0, 467.6, 468.2, 468.8, 469.3, 469.9, 470.4, 470.8, 471.2, 471.7, 472.1, 472.4, 472.8, 473.1, 473.4, 473.6, 473.9, 474.1, 474.3, 474.5, 474.6, 474.8, 474.9, 474.9, 475.0, 475.0 ] x_3 = np.arange(2012, 2100) y_3 = { 392.0, 394.0, 396.0, 398.0, 400.1, 402.3, 404.4, 406.5, 408.7, 410.9, 413.1, 415.4, 417.7, 420.0, 422.4, 424.7, 427.1, 429.5, 432.0, 434.4, 436.9, 439.5, 442.0, 444.5, 447.1, 449.8, 452.5, 455.1, 457.8, 460.5, 463.2, 466.0, 468.9, 471.7, 474.5, 477.4, 480.3, 483.3, 486.2, 489.2, 492.1, 494.9, 497.6, 500.4, 503.0, 505.6, 508.1, 510.6, 513.0, 515.4, 517.7, 519.8, 522.0, 524.1, 526.2, 528.1, 530.1, 532.0, 533.8, 535.5, 537.2, 538.8, 540.3, 541.8, 543.3, 544.7, 546.0, 547.3, 548.5, 549.7, 550.8, 551.7, 552.7, 553.6, 554.5, 555.2, 556.0, 556.7, 557.2, 557.8, 558.3, 558.7, 559.0, 559.3, 559.6, 559.8, 560.0, 560.0 } plt.figure(figsize=(12, 6)) plt.plot(x, y, 'ro', label='Data points') plt.plot(xx_extended, yy_1_extended, 'r-', label='Natural Spline') plt.plot(xx_extended, yy_2_extended, 'g-', label='Parabolic Runout Spline') plt.plot(x_1, y_1, 'b-', label='Real Data') plt.plot(x_2, y_2, 'm-', label='Future Scenario 1') plt.plot(x_3, y_3, 'c-', label='Future Scenario 2') plt.xlabel('Year') plt.ylabel('CO2 Concentration (ppm)') plt.legend(loc='best') plt.grid(True) plt.title('Global Mean CO2 Mixing Ratio over Time') plt.show()