I am trying to numerically solve some continuous-time macroeconomic models. For example, let's consider a simple growth model (as an HJB equation):
$$ \rho V(k) = \max_{c} \frac{c^{1-\sigma}}{1-\sigma} + \frac{\partial V}{\partial k}\cdot (k^\alpha-\delta k - c) $$
where $k\in(0,\infty)$ is capital, $\alpha \in(0,1)$, $\sigma>1$, $\delta\in(0,1)$ is depreciation rate, and $c>0$ is consumption of the household. Solving the RHS gives optimal consumption $c=(\frac{\partial V}{\partial k})^{-\frac{1}{\sigma}}$. Thus, the equation to solve becomes a non-linear free boundary problem (ODE in this example, but PDE in my other applications):
$$ \rho V(k) = \frac{1}{1-\sigma}(\frac{\partial V}{\partial k})^{-\frac{1-\sigma}{\sigma}} + \frac{\partial V}{\partial k}[k^\alpha - \delta k - (\frac{\partial V}{\partial k})^{-\frac{1}{\sigma}} ] $$
To well define $c$, $V$ must be strictly increasing on $k$. Meanwhile, my domain knowledge imposes a constraint that $V$ must be strictly concave. Thus, there are two shape constraints.
Also, the available boundary conditions are weak and there is singularity at $k=0$. Say: $\lim_{k\to\infty}V'(k)=0$, $\lim_{k\to0}V'(k) =\infty$ and $\lim_{k\to0}V(k)=-\infty$.
The first problem: I tried finite difference method first, and tried different feasible initial guesses (strictly increasing and strictly concave). However, the fixed point iteration often failed because, in some iteration $j$, the guess $V^j(k)$ is not strictly increasing at some grid points. My question is how to keep the shape constraints in fixed point iteration?
The second problem: When I use different initial guesses for finite difference method, even if the iteration converges, it may converge to different fixed points $V$. Does it mean that this ODE is not well-posed? If it is well-posed, could you please figure out what conditions I missed?
The third problem: I also tried other methods to solve such an ODE/PDE such as least squares. For example, suppose I have got a parametric guess $\hat{V}(k|\theta):=\theta_1 k^{\theta_2}+\theta_3$ using my domain knowledge. Then I minimize the squared HJB residual at grid points subject to the following shape constraints: $$ \begin{aligned} \text{monotonicity: }& \theta_1\theta_2 > 0 \\ \text{concavity: }& \theta_2 < 1 \end{aligned} $$ However, using different initial guesses leads to different solutions, while the minimized residuals are all very close to zero (successfully solved). It looks like the ODE is not well-posed but I could not find proper boundary conditions due to the singularity at $k=0$. Could you please suggest some possible solutions? I read some papers about dealing with singularity in linear PDEs but could not fit them into my scenario.
Thank you very much!