Write julia version of basic transition dynamics

[jlperla]

We should Julia version of one of our basic dynamics tests, such as https://github.com/econtoolkit/continuous_time_methods/blob/master/matlab/tests/discretize_time_varying_univariate_diffusion_test.m#L47

I wrote this up as https://github.com/JuliaDiffEq/PDERoadmap/blob/master/BellmanDiffusion.md which you can modify as you see fit.

• You can just replace https://github.com/JuliaDiffEq/PDERoadmap/blob/master/operator_examples/simple_stationary_HJBE_reflecting.jl and then extend it (maybe getting rid of the "stationary" part of the name?)
• Since we know how to create the basic setup of the L Q product in these simple cases, we can test this without the operator composition but still using the Julia DIffEq setup.
• Create the L Q operator as a https://docs.julialang.org/en/stable/stdlib/linalg/#Base.LinAlg.Tridiagonal
• Using the L Q at the terminal T, solve the system with u = L \ b. This will exploit the tridiagonal solver.

Given that matrix as a function of the time, we can do something along the run of

using DifferentialEquations
function f(du,u,p,t)
  L = < calculation of LQ as a tridiagonal matrix>
  A_mul_B!(du,L,u)
  du .-= c(t,L.grid)
end
solve(f, Euler(), tstops=0:0.01:10)

[jlperla]

Looking like some great progress. I am still not 100% sure about the composition of the operators.

• I thought that we would have the composition as a L_1_plus and a L_2 of size MxM instead and then multiply the constants in line with #24.
• I think it works here because of the homogeneous boundary conditions, so we can do our tweaks to the corners of the L_1 and L_2 individually and then do the multiplication and subtraction from the r I. But I could be wrong.
• It doesn't especially matter if it delivers the same L (which should be the L Q in the current notation). But I am a little worried about the r I subtraction as I am having trouble thinking through the offsets.

[jlperla]

The latest commit isn't compiling for me. the L_2 referes to the L_tilde, wich hasn't been implemented yet. Try to recompile it and you will see.

Also, can I move some of this to Julia files? Do you have atom up and running now?

[stevenzhangdx]

Sure, I haven’t upload the leatest one. The one you saw is uploaded before we fixed The issue for DifferentialEquations.jl. I have fixed it, and I start using Atom now.

On Wed, May 16, 2018 at 4:25 PM Jesse Perla notifications@github.com wrote:

The latest commit isn't compiling for me. the L_2 referes to the L_tilde, wich hasn't been implemented yet. Try to recompile it and you will see.

Also, can I move some of this to Julia files? Do you have atom up and running now?

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/JuliaDiffEq/PDERoadmap/issues/23#issuecomment-389696860, or mute the thread https://github.com/notifications/unsubscribe-auth/AdUH1_vcRovHpdRN1h8R8IRW-o_9MIeRks5tzLV-gaJpZM4T8L6F .

-- Dongxiao Zhang (Steven) Ph.D Student Vancouver School of Economics The University of British Columbia

[jlperla]

Sounds good. I pushed a variation of the code for atom that you should get as an update. You can delete it, or maybe check out the only change of substance (where I pass in the algorithm type as a setting

[jlperla]

Steven, To finish this issue, why don't you follow some of the modifications from https://github.com/JuliaDiffEq/PDERoadmap/issues/25#issuecomment-390057672

In particular, make sure to pull from the server to get the copy I put up there,

• [x] Change it to go backwards in time. It looks like this requires both swapping the order of the time range and making sure you get the sign correct on the addition/multiplication
• [x] Change the discretization of the L_1_plus and L_2 operators to avoid the resizing/etc. as we discussed. That is, directly constructing a Tridiagonal.
  • just create the 3 diagonals (of size M-1, M, and M-1 respectively) directly and then create the Tridiagonal matrix
  • I may convert these over to BandedMatrices.jl at some point later.
• [x] Test the dynamics a little bit with sufficient time variation so that we can avoid the numerical instability Chris was talking about in #25
  • In particular, try https://github.com/JuliaDiffEq/PDERoadmap/blob/master/transition_dynamics_examples/transition_dynamics_example.jl#L49 with a smaller number for the division by t so that there is large enough d/dt term in the system.
• [x] Start with adaptive time-stepping if there is enough time-variation
  • But you can try out a few of the algorithms Chris points out in his comment to get a feel for stability.
  • If things are unstable, then feel free to try a fixed dt size.

Close the issue and push to github when you think it is complete, and we can discuss next steps

[stevenzhangdx]

I changed our function to go backwards in time and tested different solving algorithm for stability. I wrote a few comments to show the stability of each algorithm. Most of them are quite stable with our numerical setup, but I could not figure out how to properly plot graphs. It seems like the plotting requests to adjust the corresponding time axis.

[jlperla]

Great! The code looks very clean. I think the plotting problems are a bug: https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/300 but could be wrong.

I will close the bug and we can do the plotting/etc. separately.