<jupyter_start><jupyter_text>Домашняя лабораторная работа №4 по вычислительной математикеДерж<NAME>, Б01-909 группа  Задача __VI.9.32__  Метод Ньютона<jupyter_code>import numpy as np
import pandas as pd
pUSA = {
    1910: 92_228_496,
    1920: 106_021_537,
    1930: 123_202_624,
    1940: 132_164_569,
    1950: 151_325_798,
    1960: 179_323_175,
    1970: 203_211_926,
    1980: 226_545_805,
    1990: 248_709_873,
    2000: 281_421_906
}

years = list(pUSA.keys())
pop = list(pUSA.values())

pop2010 = 308_745_538
def fillB(x, f):
  assert len(x) == len(f), 'x and f are not equal!!'
  b = [[0] * (len(f) - i) for i in range(len(f))]

  for i, _ in enumerate(f):
    if i == 0:
      for j in range(len(f)):
        b[0][j] = f[j]
    else:
      for j in range(len(b[i])):
        bi1 = b[i - 1]
        b[i][j] = (bi1[j + 1] - bi1[j]) / (x[j + i] - x[j])
  return b

b = fillB(years, pop)

df = pd.DataFrame(data=b, index=[f'b{i}' for i in range(len(b))])
df = df.transpose()
df['x'] = years
cols = list(df.columns)
cols = cols[10:] + cols[:10]
df = df.reindex(columns=cols)
df<jupyter_output><empty_output><jupyter_text>Так как 2010 ближе к нижнему концу таблицы, возьмем нижние числа (нижняя диагональ)<jupyter_code>n_coeffs = [bi[-1] for bi in b]
n_coeffs

def Newt(x, data_x, cfs):
  res = 0
  for i, bi in enumerate(cfs):
    mul = bi
    for j in range(i):
      mul *= x - data_x[len(data_x) - j - 1]
    res += mul
  return res

prog = int(Newt(2010, years, n_coeffs))
print(f'Прогноз = {prog}')
print(f'Delta = {prog - pop2010}')<jupyter_output>Прогноз = 827906509
Delta = 519160971<jupyter_text>Метод сплайнов Используем кубический сплайн вида $$S(x) = a_i + b_i (x- x_i) + \frac{c_i}{2} (x-x_i)^2 +  \frac{d_i}{6} (x-x_i)^3$$ Зададим граничные условия как $$f''(1910) = f''(2000) = 0$$<jupyter_code>from numpy.polynomial.polynomial import Polynomial as Poly
from matplotlib import pyplot as plt<jupyter_output><empty_output><jupyter_text>Нахождение коэффициентов сплайна осуществляется через коэффициенты $c_i$, являющимися решениями уравнения $Mx = f$, где$$M = \left(\begin{matrix}C_1 & B_1 & 0 & 0 & \dots & 0 & 0\\A_2 & C_2 & B_2 & 0 & \dots & 0 & 0\\0 & A_3 & C_3 & B_3 & \dots & 0 & 0\\\dots & \dots & \dots & \dots & \dots & \dots & \dots\\\dots & \dots & \dots & \dots & \dots & \dots & B_{n-1}\\C_1 & B_1 & 0 & 0 & \dots & A_n & C_n\\\end{matrix}\right), где \: n=8$$Граничные условия выражаются как $$c_0 = c_{10} = 0$$$$A_i = x_i - x_{i-1}, i = \overline{2\dots n}, \: \:B_i = x_{i+1} - x_{i}, i = \overline{1\dots n-1}, \:\:C_i = 2\cdot\left(x_{i+1} - x_{i-1}\right), i = \overline{1\dots n}$$$$F_i = 6 \cdot \left(\frac{y_{i+1} - y_i}{x_{i+1} - x_i} - \frac{y_{i} - y_{i-1}}{x_{i} - x_{i-1}}\right)$$<jupyter_code>N = len(years) - 2
M = np.zeros(shape=(N, N))
for i in range(N):
  if i != N - 1:
    M[i][i + 1] = years[i + 2] - years[i + 1]
  M[i][i] = 2 * (years[i + 2] - years[i])
  if i != 0:
    M[i][i - 1] = years[i + 1] - years[i]

b_col = np.zeros(shape=(N, 1))
for i, elem in enumerate(b_col):
  h_i = years[i + 1] - years[i]
  h_i1 = years[i + 2] - years[i + 1]
  b_col[i] = 6 * ((pop[i + 2] - pop[i + 1]) / h_i1 - (pop[i + 1] - pop[i]) / h_i) 
def run_through(matr, b):
  # tun though coeffs
  assert b.shape[0] == matr.shape[1], "Ivalid sizes"

  sz = b.shape[0]
  x = np.zeros(sz)
  alpha = np.zeros(sz)
  beta = np.zeros(sz)

  c = matr[0][0]
  alpha[0] = - matr[0][1] / c
  beta[0] = b[0] / c

  # direct  run through
  for k in range(1, sz):
    c = matr[k][k] + matr[k][k - 1] * alpha[k - 1]
    if k != sz - 1: 
      alpha[k] = - matr[k][k + 1] / c
    beta[k] = (b[k] - matr[k][k - 1] * beta[k - 1]) / c

  x[-1] = beta[-1]

  for k in range(sz - 2, -1, -1):
    x[k] = alpha[k] * x[k + 1] + beta[k]
  
  return x

coeffs = np.zeros(shape=(N + 1, 4))
sol = run_through(M, b_col)

for k, coef in enumerate(coeffs):
  hk = years[k + 1] - years[k]
  coef[0] = pop[k + 1]
  ck = coef[2] = sol[k] / 2 if k < len(sol) else 0
  ck_prev = coeffs[k - 1][2] if k != 0 else 0

  coeffs[k][3] = (ck - ck_prev) / (3 * hk)
  coeffs[k][1] = (pop[k + 1] - pop[k]) / hk + (2 * ck + ck_prev) * hk / 3

polys = [Poly(cfs) for cfs in coeffs]
plt.figure(figsize=[27, 8])
Y = []
X = []
for i, pol in enumerate(polys):
  z = np.linspace(-years[i + 1] + years[i], 0, num=20)
  Y += list(pol(z))
  X += list(np.linspace(years[i], years[i + 1], num=20))

plt.title('Сравнение методов интерполяции')
plt.plot(X, Y, label='Сплайн интерполяция')
plt.plot(years, pop, 'o', label='Исходные данные')
newton_y = [Newt(x, years, n_coeffs) for x in X]
plt.plot(X, newton_y, '--r', label='Метод Ньютона')
plt.legend()
plt.grid()
calc2010 = polys[-1](2010 - 2000)
print(f'Исходное значение за 2010: {pop2010}')
print(f'Посчитанное методом сплайн за 2010: {calc2010}')
print(f'Посчитанное методом Ньютона за 2010: {prog}')<jupyter_output>Исходное значение за 2010: 308745538
Посчитанное методом сплайн за 2010: 314133939.0
Посчитанное методом Ньютона за 2010: 827906509
{'fraction_non_alphanumeric': 0.19530309112806102, 'fraction_numerical': 0.0792854275391409, 'mean_word_length': 3.262617621899059, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '584490', 'license': 'MIT', 'n_tokens_mistral': 2661, 'n_tokens_neox': 2426, 'n_words': 563, 'path': 'Lab4/Lab4.ipynb', 'repo_name': 'derzhavin3016/CompMath', 'size': '4987'}