Facebook
From hi, 7 Months ago, written in Python.
Embed
Download Paste or View Raw
Hits: 139
  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.