I have a function of the following type where $p(i)$ is a variable on which the agent can optimize and $p(-i)$ a vector of variables that he has to take as parameters.

demand function

This is a demand function and from this function, one can compute a profit function.

My goal is to find the optimal solutions $p_1^*,...,p_N^*$ by computing the first-order conditions of each profit function with respect to $p(i)$ and setting each of them to 0.

Following the first answer of @user42582, I modified my code as follows:

(* Demand function *)
(* Takes a vector p1,...,pn and optional scalar arguments v and g and
returns a vector q1, ..., qn *)
(* Rn to Rn *)
q[p_, v_: V, g_: \[Gamma]] := Module[
  {meanP = Mean[p]},
  v - p - g (p - meanP)]
(* Profit function *)
(* Takes a vector p1,...,pn and optional scalar arguments cost, v and \
g and returns a vector profit1,...,profitn *)
(* Rn to Rn *)
profit[p_, cost_: c, v_: V, g_: \[Gamma]] := (p - cost) * q[p, v, g]

(* Getting the Jacobian matrix *)
(* Differentiating a Rn to Rn function wrt to X => N x N matrix *)
dvector[p_, cost_: c, v_: V, g_: \[Gamma]] := 
 D[profit[p, cost, v, g], {p}]
(* Displaying the Jacobian on an example *)
MatrixForm[dvector[{x, y, z}]]

(* Displaying only the element 2,2 of the Jacobian matrix *)
dvector[{x, y, z}][[2, 2]]

(* Solving for the prices *)
solution[p_, cost_: c, v_: V, g_: \[Gamma]] := Module[
  { n = Length[p], foc, leqs},
  (* Get the derivatives *)
  foc = Diagonal[dvector[p, cost, v, g]];
  (* Set the equations *)
  leqs = Thread[foc == ConstantArray[0, n]];
  (* Solve the equations *)
  Solve[leqs, p][[-1]]

(* Check on an example *)
solution[{x, y, z}]

Now, the solution function is fully functionnal on any example price vector of defined length like {x,y,z} or {p1, p2, p3, p4, p5} for instance.

Yet, my goal is to find the general solution of the system, which is a function of N the length of my price vector. More precisely, if p(1),...,p(N) is my price vector, I would like to express, for any i, p(i) as a function of the primitives of the model as well as N, the length of the price vector.

Is there a way to solve this in its general form ?

Thanks in advance,


  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$
    – user42582
    Feb 4, 2019 at 22:12
  • $\begingroup$ I took the Tour, thanks :)! $\endgroup$
    – MattL
    Feb 4, 2019 at 22:56
  • 1
    $\begingroup$ MattL, you can use solution in soln[n_Integer] := Module[{prices = Array[Subscript[p, #] &, n]}, solution[prices]] to get a function that takes the length of the price vector as input (and use it as soln[3], soln[6] etc.) $\endgroup$
    – kglr
    Feb 4, 2019 at 23:38
  • 1
    $\begingroup$ Youi ask for a solution as a function of "...as well as N, the length of the price vector". Why do you think this is possible? As a start, have you found solutions for $N=1$, $N=2$, or $N=3$, for example? $\endgroup$
    – Somos
    Feb 5, 2019 at 1:10
  • $\begingroup$ you don't need to modify your question unless you need to add new information $\endgroup$
    – user42582
    Feb 5, 2019 at 7:45

1 Answer 1


In Mathematica scalars eg $v$ and vectors can be mixed gracefully, courtesy of the beautiful kernel; this means, practically, that if eg $p$ is a vector or-more loosely speaking (this, too, is a distinction that is sometimes blurred but certainly not without utility)-a List, then something like v+p evaluates into the equivalent of $(v+p_1,...,v+p_n)$.

Also, when summing over a list, Total is a more suitable choice instead of Sum; the later makes use of indexes (the i part) which is not needed in order to access a list. Mathematica, unlike traditional procedural languages, does not need (it surely can, but it's not compelled) to use indexing of lists. In the present context, it wouldn't do much harm I guess but having to deal with the extra indexes makes it a bad practice to sum over immutable lists like that.

So, having said that, the following is an implementation of quantite:

  Clear[V, γ]
  q[p_, v_: V, g_: γ] := Module[{meanP = Mean[p]},
    V - p - \[Gamma] ( p - meanP )

The argument of q is intended to be a vector of prices and the output is a vector of quantities. Evaluate q[{a,b}] to verify that for the price vector {a, b} the output is as expected.

Similarly, define the profit function as

profit[p_, const_:c, v_: V, g_: γ] := (p - const) q[p, v, g]

Again, evaluate profit[{a, b}] to make sure the output is appropriate for input price vector {a, b}.

Finally, formulating the foc equations does not need to make use of indexing.

dvector[p_, const_: c, v_: V, g_: γ] := D[profit[p, const, v, g], {p}]

Admittedly, this implementation of dvector is not the most efficient; the derivatives (actually, it's the Jacobian of profit) are recalculated for every different price vector; a more robust implementation would probably first calculate the derivatives and store the result for later use; the chosen implementation is coded faster and is more transparent to someone with little experience with Mathematica.

At this point, please make sure you understand how D does indeed produce the Jacobian.

Solving for the prices is fairly straightforward:

solution[p_, const_: c, v_: V, g_: γ] := Module[{n = Length[p], foc, veqs},
  (* obtain the derivatives *)
  foc = Diagonal[dvector[p, const, v, g]];
  (* set-up the equations *)
  veqs = Thread[foc == ConstantArray[0, n]];

  Solve[veqs, p][[-1]]

Make sure to evaluate solution[{a, b}] to verify the solution is as expected.

  • $\begingroup$ Thanks a lot for this answer! I learned a lot from reading and rearranging your script. Among things I learned : optional arguments using the _: pattern, named character like [Gamma], threading list to equations like Thread[{a,b}=={x,y}] which returns {a==x, b==y}, that D can differentiate R^n to R^n functions with respect to a vector to get the Jacobian etc etc. Yet, it does not answer my initial question. How can I solve this system for its general solution in function of N, the length of my price vector ? $\endgroup$
    – MattL
    Feb 4, 2019 at 22:52
  • $\begingroup$ I don't know if my question is clear but I would like to get the general solution of the system describing the i-th element of my solution price vector as a function of the primitives of the model (const, v, g and N the length of the price vector) $\endgroup$
    – MattL
    Feb 4, 2019 at 22:54
  • $\begingroup$ it's nice to know you found something helpful; I'll be honest with you, I'm not sure if I'm answering your homework and that's something I'd hate to do; my answer irons out some rough patches in your code; take a look at the comments by kglr and Somos and consider how using induction would eventually help $\endgroup$
    – user42582
    Feb 5, 2019 at 7:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.