Лабораторна робота №2

з дисципліни «Математичне моделювання економічних та екологічних процесів»

Варіант №11

Завдання 11. 1:

Економіка країни розбита на дві виробничі галузі (промисловість та сільське господарство). За минулий рік повний випуск промислових виробництв у вартісній формі був розподілений таким чином:

1120 млн. грн. для виробничих потреб промисловості;

950 млн. грн. для виробничих потреб сільського господарства;

1970 млн. грн. для споживання населення (згідно попиту на цю продукцію).

В той же час повний випуск сільськогосподарської продукції (у вартісній формі) був розподілений таким чином:

850 млн. грн. для виробничих потреб промисловості;

720 млн. грн. для виробничих потреб сільського господарства;

1250 млн. грн. для споживання населення (згідно попиту на цю продукцію).

На наступний рік прогнозується зміна попиту населення на вітчизняну продукцію, в т. ч. на промислові вироби до 2150 млн. грн та на сільськогосподарську продукцію до 1530 млн. грн. Який повний випуск промислової продукції та повний випуск сільськогосподарської продукції зможуть задовольнити новий попит?


In [1]:
import numpy as np

In [2]:
# Дані промисловості
x11 = 1120
x12 = 950
y1 = 1970

#Дані сільського господарства
x21 = 850
x22 = 720
y2 = 1250

x1 = x11 + x12 + y1
x2 = x21 + x22 + y2

Будуємо технологічну матрицю


In [3]:
a11 = x11 / x1
a12 = x12 / x2
a21 = x21 / x1
a22 = x22 / x2

A1 = np.array([[a11, a12], [a21, a22]])
print(A1)


[[ 0.27722772  0.33687943]
 [ 0.21039604  0.25531915]]

Будуємо матрицю моделі Леонтьєва $$ \begin{cases} x_1 = a_{11}x_1 + a_{12}x_2 + y_1 \\ x_2 = a_{21}x_1 + a_{22}x_2 + y_2 \end{cases} $$


In [1]:
A2 = A1
A2[0, 0] = 1 - A2[0, 0]
A2[0, 1] = 0 - A2[0, 1]
A2[1, 0] = 0 - A2[1, 0]
A2[1, 1] = 1 - A2[1, 1]
print("A2: ")
print(A2)

# Нові об'єми
Y = np.array([2150, 1530])

X = np.linalg.solve(A2, Y)
print("X: ", X)


---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-2319d5f3e2f8> in <module>()
----> 1 A2 = A1
      2 A2[0, 0] = 1 - A2[0, 0]
      3 A2[0, 1] = 0 - A2[0, 1]
      4 A2[1, 0] = 0 - A2[1, 0]
      5 A2[1, 1] = 1 - A2[1, 1]

NameError: name 'A1' is not defined

Отже, повний випуск промислової галузі дорівнює 4528.63930886 млн. грн., а сільського господарства 3334.05615551 млн. грн.

Завдання 11. 2:

$$ A = \begin{pmatrix} 0.15 & 0.55 & 0.1 \\ 0.35 & 0.35 & 0.2 \\ 0.15 & 0.15 & 0.2 \end{pmatrix} $$

In [5]:
A = np.array([[0.15, 0.55, 0.1], [0.35, 0.35, 0.2], [0.15, 0.15, 0.2]])
print(A)


[[ 0.15  0.55  0.1 ]
 [ 0.35  0.35  0.2 ]
 [ 0.15  0.15  0.2 ]]

Знайдемо коефіцієнти характеристичного поліному


In [6]:
coefs = np.poly(A)
print(coefs)


[ 1.    -0.7   -0.085  0.016]

Знаходимо вектор v власних чисел матриці А


In [7]:
v, e = np.linalg.eig(A)
print(v)


[ 0.78249579 -0.19007257  0.10757678]

Знайдемо число Фробеніуса, яке є найбільшим дійсним додатним серед простих коренів характеристичного рівняння матриці А.


In [8]:
frob = max(v)
print(frob)


0.782495786753

Для знаходження правого вектора Фробеніуса потрібно вибрати з матриці V власних векторів той вектор, що відповідає числу Фробеніуса у матриці E.


In [9]:
left_frob = e[:, np.argmax(v)]
print(left_frob)


[-0.64680812 -0.68162645 -0.34208863]

Для знаходження лівого вектора Фробеніса потрібно повторити дії знаходження власних векторів та виведення вектора, відповідного числу Фробеніуса, але уже для транспонованої матриці до A.


In [10]:
A_trans = np.transpose(A)
v_t, e_t = np.linalg.eig(A_trans)
right_frob = e_t[:, np.argmax(v_t)]
print(right_frob)


[-0.51567456 -0.77927958 -0.3560942 ]

Знайдемо матрицю повних витрат B. $$ B = (E - A)^{-1} $$


In [11]:
B = np.linalg.inv(np.eye(A.shape[0]) - A)
print(B)


[[ 2.12121212  1.96969697  0.75757576]
 [ 1.34199134  2.87878788  0.88744589]
 [ 0.64935065  0.90909091  1.55844156]]

Для того, щоб дослідити на збіжність суму ряду $$ E + A + A^2 + ... + A^N $$ до матриці повних витрат B, потрібно організувати цикл знаходження N, у якому поступово перевіряти, чи максимальний елемент матриць не менший 0.01.


In [12]:
temp_A = B - np.ones(A.shape) - A
A_in_N = A
N = 1
while np.amax(temp_A) > 0.01:
    A_in_N = A_in_N.dot(A)
    temp_A = temp_A - A_in_N
    N = N + 1
print("N: ", N)


N:  22

Знайдемо вектор цін, якщо вектор доданої вартості в цінах s = (0.25, 0.35, 0.3). $$p = p*A + s$$ $$ p(E - A) = s $$ $$ p = s * (E - A) ^{-1} $$ $$ p = s * B $$


In [13]:
s = np.array([0.25, 0.35, 0.3])
p = s.dot(B)
print(p)


[ 1.19480519  1.77272727  0.96753247]