Least-squares, Fourier series

Least-squares

As far as I know, Sage does not have a built-in method to find a “least-squares solution” to a system of linear equations. The description of a least squares solution to Ax=b as a solution to A^TAx=A^Tb is easy to work with in Sage. For example, to find a least-squares solution to the system

x₁ + x₂ = 3
2x₁ + 3x₂ = 1
x₁ – x₂ = 5

we can type:
A = matrix([[1, 1],[2,3],[1,-1]])
B = A.transpose()
v = vector([3,1,5])
(B*A).solve_right(B*v)

or, more succinctly,
A = matrix([[1, 1],[2,3],[1,-1]])
(A.transpose()*A).solve_right(A.transpose()*vector([3,1,5]))

We could also write a function that finds least-squares solutions for us:

def least_squares(A,v):
     return (A.transpose()*A).solve_right(A.transpose()*v)

(Again, indentation is important: this is Python.) Then we can call it with the matrix and vector of your choice:
B = matrix([[1,2],[3,4],[5,6]])
v = vector([7,5,9])
B.solve_right(v)
least_squares(B,v)

(There is also built-in code in Sage and numpy for fitting curves to data.)

Fourier series

Sage has some rudimentary support for Fourier series, as part of the “piecewise-defined function” class, but it seems to be very slow and very flaky. Let’s implement our own. To make things run reasonably efficiently, we’re going to have Sage do numerical, rather than symbolic, integrals.

As a warm up, to integrate x^2 from 0 to 1, you type:
numerical_integral(x^2,0,1)
The first entry is the answer, while the second is an error estimate.

With this in hand, here’s a function to return a Fourier approximation:

def numerical_fourier(f, n):
    ans = 0
    for i in range(1,n+1):
        ans+= (numerical_integral(f(x)*sin(i*x), -pi,pi)[0])*sin(i*x)/float(pi)
        ans+= (numerical_integral(f(x)*cos(i*x), -pi,pi)[0])*cos(i*x)/float(pi)
    ans+= numerical_integral(f(x),-pi, pi)[0]/(2*float(pi))
    return ans

The inputs to this function are the function you want to approximate, and the number n of sine terms to use.

Let’s approximate the function f(x)=e^sin(x), which is a random-ish periodic function of period 2pi, by a Fourier series with 7 terms:
numerical_fourier(lambda x:exp(sin(x)), 3)


This “lambda x: exp(sin(x))” is a way of quickly defining and passing a function, called lambda calculus.

The answer is not so enlightening. It’s a bit better if we plot it:
plot(exp(sin(x)),(x,-2*pi,2*pi), color='red')
plot(numerical_fourier(lambda x:exp(sin(x)), 3), (x,-2*pi,2*pi))


Let’s plot several Fourier approximations on the same plot, to watch them converging:
plot(exp(sin(x)),(x,-2*pi,2*pi), color="red")+plot(numerical_fourier(lambda x:exp(sin(x)), 1), (x,-2*pi,2*pi), color='green')+plot(numerical_fourier(lambda x:exp(sin(x)), 2), (x,-2*pi,2*pi), color='blue')+plot(numerical_fourier(lambda x:exp(sin(x)), 3), (x,-2*pi,2*pi), color='orange')

It converges pretty quickly — it’s hard to even see the difference between orange (approximation) and red (exact).