Solve variation with non-uniform grid

[jlperla]

Base this off of #8, though it could easily be done with the LCP.m based one.

Essentiallly, just do the finite differences with a fixed, but non-uniform grid. The key step will be that dx is a vector of first differences of the x vector (provided by the user).

Then when creating the X,Y,Z matrices or generating the $A$ matrix, care will need to be taken to divide by ./dx The 2nd difference operator changes as well: Care will be need to be taken if diving by dx(i) vs. dx(i+1) vs. dx(i-1) given the current setup.

See http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf page 22 for most of the work done and for verification.

[jlperla]

Start here with the operator_discretization_finite_differences.tex modifications. We can get the matlab code written after the algebra is checked.

The complexity here is that what $\Delta$ we are dividing by changes depending on what the upwind direction is. Consider doing this in two stages: first, figure out the correct algebra without worrying about how it is structured in the matrix. This is essentially replicating (13), (14), (17) and (20). (don't worry about the non-reflecting boundary value). For this, you can leave the $\Delta$ differences in the denominator. After you get that, try to figure out how to rearrange it so we don't need to divide by the equivalent to the $\Delta^2$ term, as we did when we rearranged the uniform grid example.

I put a placeholder for where to put the text for this in section 1.3. The (29) to (31) guesses are probably wrong, and will depend on how the multiplication by the approxiate Delta vectors is done. Note my suggestion in (34) that there may be a way to come up with a $\Delta$ vector which reflects the appropriate direction from the upwind process.

[stevenzhangdx]

Here is the attempt for this non-uniform grid issue.

Based on the recommendation of second derivative choice from http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf, a good candidate approximation is \begin{equation} vi^{''} = \frac{\Delta{-}v{i+1} - (\Delta{-}+\Delta_{+})vi+\Delta{+}v{i-1}}{\frac{1}{2}(\Delta{+}+\Delta{-})\Delta{-}\Delta{+}} \nonumber \end{equation}. Then we can define $X$, $Y$ and $Z$ as \begin{equation} X = -\frac{\mu^{-}}{\Delta{-}} + \frac{\sigma^2}{2}\frac{\frac{\Delta{+}}{\frac{1}{2}(\Delta{+}+\Delta{-})}}{\Delta{+}\Delta{-}} \nonumber \end{equation} \begin{equation} Y = -\frac{\mu^{+}}{\Delta{+}}+\frac{\mu^{-}}{\Delta{-}} -\frac{\sigma^2}{\Delta{+}\Delta{-}} \nonumber \end{equation} \begin{equation} Z = \frac{\mu^{+}}{\Delta{+}} + \frac{\sigma^2}{2}\frac{\frac{\Delta{-}}{\frac{1}{2}(\Delta{+}+\Delta{-})}}{\Delta{+}\Delta{-}} \nonumber \end{equation} For the stability of numerical calculation, we can multiply every term in $A$ by $\Delta{+}\Delta{-}$ to obtain $\Delta{+}\Delta{-}A$. So that we can redefine $X, Y $and $Z$ as \begin{equation} X = -\mu^{-}\Delta{+} + \frac{\sigma^2}{2}\frac{\Delta{+}}{\frac{1}{2}(\Delta{+}+\Delta{-})} \nonumber \end{equation} \begin{equation} Y = -\mu^{+}\Delta{-}+\mu^{-}\Delta{+} -\sigma^2 \nonumber \end{equation} \begin{equation} Z =\mu^{+}\Delta{-} + \frac{\sigma^2}{2}\frac{\Delta{-}}{\frac{1}{2}(\Delta{+}+\Delta_{-})} \nonumber \end{equation}

If this attempt is correct, care will need to be taken to determine the result from discretize_univariate_diffusion.m, which is uniform grid or non-uniform grid whenever we are using the function discretize_univariate_diffusion.m.

[jlperla]

Great. Edit in the tex file and I will take a look at the algebra in detail

On Sat, Oct 21, 2017, 6:41 PM stevenzhangdx notifications@github.com wrote:

Here is the attempt for this non-uniform grid.

Based on the recommendation of second derivative choice from http://www.princeton.edu/~moll/HACTproject/HACT_Numerical_Appendix.pdf, a good candidate approximation is \begin{equation} vi^{''} = \frac{\Delta{-}v{i+1} - (\Delta{-}+\Delta_{+})vi+\Delta{+}v{i-1}}{\frac{1}{2}(\Delta{+}+\Delta{-})\Delta{-}\Delta{+}} \nonumber \end{equation}. Then we can define $X$, $Y$ and $Z$ as \begin{equation} X = -\frac{\mu^{-}}{\Delta{-}} + \frac{\sigma^2}{2}\frac{\frac{\Delta{+}}{\frac{1}{2}(\Delta{+}+\Delta{-})}}{\Delta{+}\Delta{-}} \nonumber \end{equation} \begin{equation} Y = -\frac{\mu^{+}}{\Delta{+}}+\frac{\mu^{-}}{\Delta{-}} -\frac{\sigma^2}{\Delta{+}\Delta{-}} \nonumber \end{equation} \begin{equation} Z = \frac{\mu^{+}}{\Delta{+}} + \frac{\sigma^2}{2}\frac{\frac{\Delta{-}}{\frac{1}{2}(\Delta{+}+\Delta{-})}}{\Delta{+}\Delta{-}} \nonumber \end{equation} For the stability of numerical calculation, we can multiply every term in $A$ by $\Delta{+}\Delta{-}$ to obtain $\Delta{+}\Delta{-}A$. So that we can redefine $X, Y $and $Z$ as \begin{equation} X = -\mu^{-}\Delta{+} + \frac{\sigma^2}{2}\frac{\Delta{+}}{\frac{1}{2}(\Delta{+}+\Delta{-})} \nonumber \end{equation} \begin{equation} Y = -\mu^{+}\Delta{-}+\mu^{-}\Delta{+} -\sigma^2 \nonumber \end{equation} \begin{equation} Z =\mu^{+}\Delta{-} + \frac{\sigma^2}{2}\frac{\Delta{-}}{\frac{1}{2}(\Delta{+}+\Delta_{-})} \nonumber \end{equation}

If this attempt is correct, care will need to be taken to determine the result from discretize_univariate_diffusion.m, which is uniform grid or non-uniform grid whenever we are using the function discretize_univariate_diffusion.m.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/econtoolkit/continuous_time_methods/issues/9#issuecomment-338444386, or mute the thread https://github.com/notifications/unsubscribe-auth/AMf-IbHyppBz6eRA0qSAjL8qtYX7TRqSks5sup1BgaJpZM4PSNmC .

[stevenzhangdx]

Explanation of this attempt has been added to the tex file. Please take a look at it.

[jlperla]

Looks good. A few question:

• To make more stable, you multiplied by $\Delta{+} \Delta{-}$. Was this suggested in the reference, or was it your approach?
• Could we just multiply by one of them, e.g. $\Delta{-}$, and then we would have a couple of places with $\Delta{+}/\Delta{-}$ but otherwise the whole expression would not be multiplied by any $\Delta$ in the numerator? Since $\Delta{+}/\Delta_{-} \approx 1$, I think this is likely stable.
• You mention that we need to deal with the $\Delta$ to the ghost nodes. I think your approach is perfectly reasonable.
• The last thing we need to do here is make sure that the boundary values are taken into account, so that we can have the full solution as a system.

[stevenzhangdx]

Multiplying everything by $\Delta{+}\Delta{-}$ is my approach. In terms of stability, I think your suggestion is better. I well revise the code later today and write corresponding equations into tex file.

[jlperla]

Sounds great. My goal is to treat the backwards drift as the most likely scenario in this setup, so keep that in mind when deciding if multiplying by $\Delta{+}$ or $\Delta{-}$ is more natural. For example, in (28) to (30) it looks to me like the $\mu^+$ terms would knock out most of the $\Delta{+}$, if the drift is negative, which makes me think we should multiply by $\Delta{-}$.

Also add in the boundary values for the $A$ matrix (remembering that we only need to worry about reflecting boundaries for now) and we can write the equivalent to (12) and (45). Probably something like $\Delta_{-} \mathcal{A}(v({x_i})) \approx A v$??? where that is a pointwise multiplication on the left?

[stevenzhangdx]

I updated the multiplication and finalized corresponding equations to the tex file. Please take a look at it and let me know if there is any mistake.

[jlperla]

Thanks, looking good. A few points:

• The equation (40) to (47) needs to have the general $\Delta$ rather than the particular oneone. I made a few changes, which I think are worth verifying
• In particular, check that (42) to (44) is correct given the assumption on ghost node distances $\Delta{1,-} = \Delta{1,+}$. This basically reduces to the constant $\Delta$ scenario as far as I can tell.
• Similarly, (45) to (47) simplify in a similar way and should be checked.
• Check that I did (34) to (39) correctly when adjusting for the $\Delta_{-}$ division. Need to be careful as that is a vector, so I defined a diagonalized version of this.
• Then see if (68) and (69) are correct for solving the linear HJBE given the diagonalized $\Delta_{-}$.

[jlperla]

I just added another equation which may offset the equation numbers in previous comment by 1.

One question when solving the HJBE is which $\Delta$ to use for discretized $u(x)$ and $\rho v(x)$. Try to see in the Moll code if there is a hint which one, or if ti chagnes depending on the wind.

[stevenzhangdx]

I checked modifications you made and most of them look good to me. For (67), in the non-uniform grid case, degrees do not match when we are multiplying $\Delta$ on the last row. Since $\Delta$ is just a constant multiplier on the last row, my guess is that solving this system of equations without this multiplication on last row does not affect the result. Also, for section 1.6, I could not find any suggestion for choice of $\Delta$ from Moll's code. To me, $u(x)$ and $\rho v(x)$ are discretized in terms of vector x. Which multiplier we choose does not affect the discretization of $u(x)$ and $\rho v(x)$.

[jlperla]

Yes, we might need to rethink exactly how (67) looks if we were to do a joint solution with the KFE. The reality is that the last row should be using a particular quadrature rule for non-uniform grids, so the $\Delta \omega$ might become trickier. At this point, not really necessary. though. I am closing this task, and will add in a new task for writing the code.