1
$\begingroup$

I'm trying to solve an ODE where there is an optimal control problem involving a tradeoff between the flow cost c and the drift rate g as well as value matching and smooth pasting conditions characterizing the boundary value xbar. Something like the following:

0 == c[a, x] + f[a, x]*V[x] + g[a, x]*V'[x] + h[x]*V''[x]
V[xbar] == psi[xbar]
V'[xbar] == psi'[xbar]

This is a class of problem that often appears in financial economics. I tried to implement this in Mathematica using FindMaximum and NDSolve, but encountered issues. Minimal Example below:

ClearAll["Global`*"]
c[a_, x_] = (1 - a^2/2)*x;
f[x_]     = 1;
g[a_, x_] = a;
h[x_]     = 0.5*x^2;
psi[x_]   = x - 1;
astar[X_?NumericQ, dV_?NumericQ] := 
   With[{xx = X, dvv = dV}, 
    FindMaximum[{c[a, xx] + g[a, xx]*dvv, a >= 0}, {a, 
       1/2}][[2]]]; // Quiet
ODE = {0 == c[a, x] + f[x]*V[x] + g[a, x]*V'[x] + h[x]*V''[x], 
   V[xbar] == psi[xbar], 
   V'[xbar] == psi'[xbar]} /. {a -> astar[x, V'[x]]}

I am trying to use a shooting method to evaluate whether my guess of is too high or too low relative to the condition V(0) = 0, and will eventually do this properly using FindRoot.

Trying with a guess of 10 (as in the example) results in errors. What is the appropriate way to implement this in Mathematica?

EDIT: To provide some more structure to the problem, x is assumed to take values $x \in \left[0, \overline{x}\right]$, and $a$ is restricted to positive values. In the MWE, I provided a convex cost function $c$ that has support over all nonnegative numbers, but it is equally interesting to consider something of the form $\left(a-\ln\left(1-a\right)\right)x$, which would only take values on the interval $\left[0, 1\right)$.

$\endgroup$
  • $\begingroup$ You use a definition f[x_]=1 with one argument, and in the equation you use this function with two arguments f[a,x]. How do you solve equation? $\endgroup$ – Alex Trounev Jul 25 '18 at 5:12
  • $\begingroup$ @AlexTrounev ammended to reflect the fact that f is a trivial function of only x. The error was mine. $\endgroup$ – Shffl Jul 25 '18 at 6:28
  • $\begingroup$ It is necessary to formulate the problem correctly or to indicate an article where this task is formulated. We do not even know the interval for x and a. $\endgroup$ – Alex Trounev Jul 25 '18 at 7:16
  • $\begingroup$ Based on @AlexTrounev's request, I've edited the original post to provide more structure. $\endgroup$ – Shffl Jul 25 '18 at 18:45
2
$\begingroup$

Here is the version of the running code as I understand this task

ClearAll["Global`*"]
x0 = .1; xbar = .9;
c[a_, x_] := (1 - a^2/2)*x;
f[x_] := 1;
g[a_, x_] := a;
h[x_] := 0.5*x^2;
psi[x_] := x - 1;
astar[X_?NumericQ, dV_?NumericQ] := 
   With[{xx = X, dvv = dV}, 
    FindMaximum[{c[a, xx] + g[a, xx]*dvv, a >= 0}, {a, 
       1/2}][[1]]]; // Quiet
ODE = {c[a, x] + f[x]*V[x] + g[a, x]*V'[x] + h[x]*V''[x] == 0, 
    V[xbar] == psi[xbar], 
    V'[xbar] == psi'[xbar]} /. {a -> astar[x, V'[x]]};
sol = NDSolve[ODE, V, {x, x0, xbar}];

{Plot3D[astar[x, y], {x, 0, 1}, {y, 0, 1}, Mesh -> None, 
 AxesLabel -> {"x", "y", ""}],Plot[Evaluate[V[x] /. sol], {x, x0, xbar}, AxesLabel -> {"x", "V(x)"}]}

fig1

| improve this answer | |
$\endgroup$
  • $\begingroup$ This code works, and in particular works even when I allow xbar to be a parameter rather than a predefined constant. Is there a better way of implementing the Shooting Method to impose the condition V[0] == 0 inside of NDSolve? My prior experience was to do it myself using FindRoot $\endgroup$ – Shffl Jul 25 '18 at 18:53
  • $\begingroup$ Looking at the code again, is there a reason you have a->astar, where you have changed the definition of astar to be the maximized reward, rather than the optimal control value a? $\endgroup$ – Shffl Jul 25 '18 at 20:12
  • $\begingroup$ To avoid singularity, we must exclude the point x = 0. Note that under the existing conditions and data the set FindRoot[][[2]] is empty, so I took the set FindRoot[][[1]]. We can not put three boundary conditions for one equation of the second order. There must be a free parameter to execute the third boundary condition, but it is not among the data. $\endgroup$ – Alex Trounev Jul 26 '18 at 4:42
  • $\begingroup$ Assuming I understood you correctly, the free parameter is the choice of xbar. The point regarding the set notation with FindRoot is interesting and one I wasn't aware of. $\endgroup$ – Shffl Jul 27 '18 at 17:36
  • $\begingroup$ You can enter a free parameter in any function, for example, psi[x_]:= x + k, then for k = 2.39 we get V[0]=0 $\endgroup$ – Alex Trounev Jul 27 '18 at 18:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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