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:

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)$.

  • $\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$ Jul 25, 2018 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, 2018 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$ Jul 25, 2018 at 7:16
  • $\begingroup$ Based on @AlexTrounev's request, I've edited the original post to provide more structure. $\endgroup$
    – Shffl
    Jul 25, 2018 at 18:45

1 Answer 1


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

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)"}]}


  • $\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, 2018 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, 2018 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$ Jul 26, 2018 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, 2018 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$ Jul 27, 2018 at 18:17

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.