# 2nd Order ODE with Optimal Control

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

• 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? Jul 25, 2018 at 5:12
• @AlexTrounev ammended to reflect the fact that f is a trivial function of only x. The error was mine. Jul 25, 2018 at 6:28
• 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. Jul 25, 2018 at 7:16
• Based on @AlexTrounev's request, I've edited the original post to provide more structure. Jul 25, 2018 at 18:45

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


• 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 Jul 25, 2018 at 18:53
• 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? Jul 25, 2018 at 20:12
• 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. Jul 26, 2018 at 4:42
• 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. Jul 27, 2018 at 17:36
• 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 Jul 27, 2018 at 18:17