CODEWITHSUNDEEP
pythonsympy
Davoud Taghawi-Nejad
Davoud Taghawi-Nejad

Reputation: 16776

sympy element wise division of vector

I have a price vector and a quantity vector. I would like to symbolically divide the two vectors in sympy. I can't even find a vector facility. So I am totally lost with sympy on how to create vector symbols and how to manipulate them.

Eventually what I would like use the lagrangian method to maximize is prod((x/p) **b) subjected to the constrained sum(x) = 1. I can do this with scalars, but not with vectors:

import sympy as sp
import sympy
from  sympy.abc  import,p1, p2, l, x1, x2,  b1, b2
sp.init_printing()
U = ((x1/p1)**b1)*((x2/p2)**b2)
L = U - l*(x1 + x2 - 1)
dL_dy = sp.diff(L, x1)
dL_dx = sp.diff(L, x2)
dL_dl = sp.diff(L, l)
sp.solve([dL_dy, dL_dx, dL_dl], (x1, x2, l))

Upvotes: 0

Views: 1643

Answers (2)

smichr
smichr

Reputation: 18939

IndexedBase is a nice object to use when you need a vector-like quantity:

>>> from sympy import *
>>> quant=IndexedBase('q')
>>> price=IndexedBase('p')
>>> i = var('i',integer=True)
>>> product(price[i]/quant[i],(i,1,2))
p[1]*p[2]/(q[1]*q[2])

You can differentiate wrt p[1], q[1], etc..., too. Would that help you formulate the problem more naturally?

Upvotes: 0

Stelios
Stelios

Reputation: 5521

Here is one way to do this.

import sympy as sp

Define (vector) variables and parameters:

# vector size (integer, user input):
n = 2
# vector variables:
x = sp.var('x0:'+str(n), positive = True)
y = sp.var('y0:'+str(n), positive = True)
# vector parameters:
p = sp.var('p0:'+str(n), positive = True)
q = sp.var('q0:'+str(n), positive = True)
# scalar parameters
b = sp.var('b', real = True)
c = sp.var('c', real = True)
# Lagrange multiplier for sum constraint:
l = sp.var('lambda')

Objective function: 

U = reduce(lambda xi, xj: xi * xj, [(xi/pi)**b * (yi/qi)**c for xi,pi,yi,qi in zip(x,p,y,q)],1)
U

(x0/p0)**b*(x1/p1)**b*(y0/q0)**c*(y1/q1)**c

Lagrangian:

L = U + l * (sum(x+y)-1)

KKT conditions (each list element must be equal to zero):

KKT = sp.simplify([sp.numer(sp.together(sp.diff(L, xi))) for xi in x]+\
        [sp.numer(sp.together(sp.diff(L, xi))) for yi in y] + [sp.diff(L, l)])

I have considered only the numerator of the derivatives in order to help the solver. This means that some solutions based on this approach may be invalid due to a corresponding zero denominator (they have to be manually checked).

The solution can be obtained now as 

sp.solve(KKT,sp.flatten([x,y,l]))

It appears that for general values of the parameters b and c, Sympy is unable to give a solution. However, solutions can be obtained for certain choices of these parameters. For example, for b=2 and c=2, the solution given is

[{lambda: y0**2*y1**2*(y0 + y1 - 1)**3/(4*p0**2*p1**2*q0**2*q1**2), 
  x0: -y0/2 - y1/2 + 1/2, 
  x1: -y0/2 - y1/2 + 1/2}]

Upvotes: 2

Related Questions