Logarithmic step size

[crisla]

I've been trying to replicate the original Adamopoulos and Restucia paper, and they say they use a logarithmic grid for farmers' ability with support close to zero to 100.

How do you change the formula for the step size in the solvers' code so that the step size is logarithmic instead of uniform (ie (xhigh - xlow) / (knots-1))? I've been searching for a while with little success so far.

[davidrpugh]

@crisla

Sounds like what you need is the ability to pass in an array of points at which you would like the solver to return the solution. To get the grid you want you can probably want to use numpy.logspace. Something like:

>>> import numpy as np
>>> xlow = 1e-3
>>> xhigh = 1e2
>>> start = np.log10(xhigh)
>>> stop = np.log10(xlow)
>>> grid = np.logspace(start, stop, 50, base=10, endpoint=True)
>>> grid
array([  1.00000000e+02,   7.90604321e+01,   6.25055193e+01,
         4.94171336e+01,   3.90693994e+01,   3.08884360e+01,
         2.44205309e+01,   1.93069773e+01,   1.52641797e+01,
         1.20679264e+01,   9.54095476e+00,   7.54312006e+00,
         5.96362332e+00,   4.71486636e+00,   3.72759372e+00,
         2.94705170e+00,   2.32995181e+00,   1.84206997e+00,
         1.45634848e+00,   1.15139540e+00,   9.10298178e-01,
         7.19685673e-01,   5.68986603e-01,   4.49843267e-01,
         3.55648031e-01,   2.81176870e-01,   2.22299648e-01,
         1.75751062e-01,   1.38949549e-01,   1.09854114e-01,
         8.68511374e-02,   6.86648845e-02,   5.42867544e-02,
         4.29193426e-02,   3.39322177e-02,   2.68269580e-02,
         2.12095089e-02,   1.67683294e-02,   1.32571137e-02,
         1.04811313e-02,   8.28642773e-03,   6.55128557e-03,
         5.17947468e-03,   4.09491506e-03,   3.23745754e-03,
         2.55954792e-03,   2.02358965e-03,   1.59985872e-03,
         1.26485522e-03,   1.00000000e-03])
>>> 

Once we pass this grid into the solver, then we need to add a method for computing desired step size which gets called at each step of the integration. Probably something like the following would work...

>>> np.diff(grid, 1)
array([ -2.09395679e+01,  -1.65549129e+01,  -1.30883856e+01,
        -1.03477342e+01,  -8.18096341e+00,  -6.46790502e+00,
        -5.11355366e+00,  -4.04279762e+00,  -3.19625327e+00,
        -2.52697164e+00,  -1.99783470e+00,  -1.57949675e+00,
        -1.24875695e+00,  -9.87272643e-01,  -7.80542018e-01,
        -6.17099892e-01,  -4.87881841e-01,  -3.85721492e-01,
        -3.04953078e-01,  -2.41097221e-01,  -1.90612505e-01,
        -1.50699070e-01,  -1.19143336e-01,  -9.41952363e-02,
        -7.44711608e-02,  -5.88772215e-02,  -4.65485858e-02,
        -3.68015130e-02,  -2.90954352e-02,  -2.30029768e-02,
        -1.81862529e-02,  -1.43781301e-02,  -1.13674118e-02,
        -8.98712488e-03,  -7.10525977e-03,  -5.61744907e-03,
        -4.44117951e-03,  -3.51121571e-03,  -2.77598231e-03,
        -2.19470361e-03,  -1.73514216e-03,  -1.37181089e-03,
        -1.08455962e-03,  -8.57457520e-04,  -6.77909620e-04,
        -5.35958275e-04,  -4.23730928e-04,  -3.35003503e-04,
        -2.64855217e-04])

...which means that the new method would look something like...

def _compute_step_size(self, grid):
    return np.diff(grid, 1)

This method should then get called outside the while loop to compute the step sizes.

step_sizes = self._compute_step_size(grid)

Then within the while loop one would just select the next step size in the sequence where appropriate.

[crisla]

That sounds great! Do you think this maybe has to do with the fact that the collocation solver was not so good at points close to zero? Maybe we needed to feed it more small values?

[davidrpugh]

@crisla

I don't think this is the problem with the collocation solver. For the most part I was using Chebyshev polynomials which automatically leads to more knots being placed near the "edges" of the domain. I will get back to the collocation solver later this week.

Let me know how you get on with using the new functionality I added in PR #21.