[> restart; with(LinearAlgebra):
[> p[1, 0]: =4: p[2, 0]: =5: M_[0]: =400:
[> Vector_p: =Matrix(1, 2, [[p[1], p[2]]]):
[> U: =3*ln(Q[1]-10)+4*ln(Q[2]-20);
/находим предельные полезности, градиент функции полезности, матрицу Гессе (матрицу вторых частных производных). Вычисляем обратную матрицу к матрице Гессе/
[> MU[1]: =simplify(diff(U, Q[1])): MU[2]: =simplify(diff(U, Q[2])):
Gradient_U: =[MU[1], MU[2]];
[> G[1, 1]: =simplify(diff(MU[1], Q[1])):
G[1, 2]: =simplify(diff(MU[1], Q[2])):
G[2, 1]: =simplify(diff(MU[2], Q[1])):
G[2, 2]: =simplify(diff(MU[2], Q[2])):
[> Gesse: =Matrix(2, 2, [[G[1, 1], G[1, 2]], [G[2, 1], G[2, 2]]]);
Gessian: =Determinant(Gesse);
Gesse_Inverse: =MatrixInverse(Gesse);
[> p[1]: =p[1, 0]: p[2]: =p[2, 0]: M: =M_[0]:
[> g: =p[1]*Q[1]+p[2]*Q[2]-M; L: =U-lambda*g;
[> sys: ={diff(L, Q[1])=0, diff(L, Q[2])=0, diff(L, lambda)=0}:
[> Local_Ravnoves: =fsolve(sys, {Q[1], Q[2], lambda});
[> lambda: =rhs(Local_Ravnoves[1]);
Q[1]: =rhs(Local_Ravnoves[2]); Q[2]: =rhs(Local_Ravnoves[3]);
/вычисление вектора эффекта дохода, вычисляем коэффициент /
[> p[1]: =p[1, 0]: p[2]: =p[2, 0]: M: =M_[0]:
[> b: =MatrixMatrixMultiply(Gesse_Inverse, Transpose(Vector_p)):
mu: =-1/MatrixMatrixMultiply(Vector_p, b);
[> C1: =MatrixMatrixMultiply(Gesse_Inverse, Transpose(Vector_p)):
[> Vector_Diff_Q_M: =MatrixScalarMultiply(C1, (-1)*mu[1, 1]);
[> Vector_Diff_Q_M_Q1: =MatrixScalarMultiply(Vector_Diff_Q_M, Q[1]);
Vector_Diff_Q_M_Q2: =MatrixScalarMultiply(Vector_Diff_Q_M, Q[2]);
/вычисляем вектор эффекта компенсации (по расчетной формуле (7.8))/
[> D1: =MatrixMatrixMultiply(Vector_p, Gesse_Inverse):
[> D2: =MatrixMatrixMultiply(C1, D1):
[> D3: =MatrixScalarMultiply(D2, mu[1, 1]):
[> D4: =MatrixAdd(D3, Gesse_Inverse):
[> D5: =Matrix(2, 1, [D4[1, 1], D4[2, 1]]):
[> D6: =Matrix(2, 1, [D4[1, 2], D4[2, 2]]):
[> Vector_Diff_Q_p1_Comp: =MatrixScalarMultiply(D5, lambda);
Vector_Diff_Q_p2_Comp: =MatrixScalarMultiply(D6, lambda);
|