Week July 19 - 23

This week I have learned how projections work in detail and wrote it up here:

http://theoretical-physics.net/dev/src/math/la.html

including all proofs (that orthogonal projection finds the closest vector and so on). At the end of the notes I have calculated some examples in 1D, so that one can see that indeed it doesn't depend on the basis and that the basis doesn't even have to be orthogonal. Then one has to use this approach to calculate the coefficients:

http://theoretical-physics.net/dev/src/math/la.html#nonorthogonal-basis


I have then implemented it in h1d. Due to the lack of time, I am now developing everything in my private branch. I'll obtain the permission in about 2 or 3 weeks, so then I'll push it into the master.

Besides that I have implemented Chebyshev points for orders greater than 48, for which I don't have exact Fekete points anymore (it'd be just a matter of running my sympy script longer, but I was hitting some accuracy issues when solving those large polynomials numerically -- one needs to obtain all roots, so Chebyshev points are ok for now). So I can now represent arbitrary polynomials in 1D.

I have implemented powering of the discrete function, it automatically determines which polynomial order it has to use and creates a new discrete function (the power) on the new mesh. I wrote lots of tests for that, and I hit an interesting bug, that my naive comparison code:

assert abs(x-y) < eps

was not good enough anymore for larger numbers and I had to read some documentation, and implement the following function:

@vectorize(0, 1)
def feq(a, b, max_relative_error=1e-12, max_absolute_error=1e-12):
    a = float(a)
    b = float(b)
    # if the numbers are close enough (absolutely), then they are equal
    if abs(a-b) < max_absolute_error:
        return True
    # if not, they can still be equal if their relative error is small
    if abs(b) > abs(a):
        relative_error = abs((a-b)/b)
    else:
        relative_error = abs((a-b)/a)
    return relative_error <= max_relative_error

Then I implemented the global H1 and L2 projections, so far the projected function is hardwired, I still need to allow the user to specify any discrete function to be projected. I need to precalculate it and so on.

I wrote bunch of tests for the projections and powers and I always discovered some bugs by writing more tests, so the progress is slow, but at least I can trust the code that is tested.

I also helped Pavel to fix couple segfaults, as well as some other things.