Polynomial Approximations with Gram-Schmidt
A friend recommended Linear Algebra Done Right, which has within it an exercise in approximating sin(x) using orthogonal projections onto a P(5) basis generated using the Gram-Schmidt procedure.
I was curious about this approximation, and how closely the derivative of the approximation would approach cos(x). 21 pages of notes and an excel spreadsheet later, I have some observations (and shortcuts) to share.
In the text, the problem is to approximate sin(x) with a 5th order polynomial over the range [-π,π].
The first step is to invoke Gram-Schmidt to find a othogonal basis to work in. The procedure is designed to produce an orthonormal basis, but I quickly found that the various square roots lying about weren't worth the clutter.
<x^n,e[j]>e[j]
= (B[j]^-2)<x^n,B[j]e[j]>B[j]e[j]
= (B[j]^-2)<x^n,v[j]>v[j]
I never seemed to need B[j] alone; it is always squared, so I just kept it that way as I worked.
Also, while the text starts the basis from e[1], I found it less confusing to start with e[0] - the subscript matching the highest order in X.
We're trying to find a basis to span P5, so it makes sense to start with the straightforward (1,x,x^2,x^3,x^4,x^5). 1 is even, and choosing that for the first vector is going to produce a basis which alternates between even and odd functions.
Since the inner product must be zero, it follows that the integrals of each v[j] must be identically zero. This isn't useful for discovery, really, but makes for a good check that you haven't mangled the symbols during the manipulation.
Similarly, as we begin calculating, we discover there are lots of π's about. The check here is that each term in v[j] must be proportional to x^mπ^n, where j = m+n. This means that, in practice, you can throw the π's away when calculating, and stick them in at the end.
What finally gave me control of discovery was to precalculate, for each v[j], the constant that sits out in front of it.
(B[j]^-2)
= A[j](n)
Again, the π's will take care of themselves, so you can instead write
B[j]^2 = π^(2j+1) B'[j]^2
A[j](n) = π^(n-j) A'[j]
but I admit that in my notes I kept the π's around as a security blanket.
With these tools at the ready, the caculation is fairly straightforward
v[n] = x^n - sum( A[j](n)v[j] )
where j ranges over each of the integers smaller than n with the same parity.
Now you can, if you are so inclined, simplify at this point, expressing v[n] as a conventional polynomial. For example, if n = 4
v[4]
= x^4 - sum( A[j](4)v[j] )
= x^4 - A[2](4)v[2] - A[0](4)v[0]
= x^4 - A[2](4)(x^2-A[0](2)v[0]) - A[0](4)v[0]
= x^4 - A[2](4)x^2 + (A[2](4)A[0](2) - A[0](4))v[0]
= x^4 - A[2](4)x^2 + (A[2](4)A[0](2) - A[0](4))
Perfectly valid, but not convenient to work with when you try to determine the norm.
Of course, I later discovered that there's a simply trick to finding the right B[j] values. Consider this expression: x^n - A[j](n)v[j] - A[k](n)v[k]
To find the norm, we square it, then take the integral. The square looks like
x^2n
-2 x^nA[j](n)v[j] + A[j](n)^2v[j]v[j]
-2 x^nA[k](n)v[k] + A[j](n)^2v[k]v[k]
+2 A[j](n)A[k](n)v[j]v[k]
and when we integrate, we get
(x^(2n+1))/(2n+1)
- 2A[j](n)<x^n,v[j]> + A[j](n)^2<v[j],v[j]>
- 2A[k](n)<x^n,v[k]> + A[j](n)^2<v[k],v[k]>
+2 A[j]A[k]<v[j],v[k]>
We've already evaluated those inner products. The last term vanishes, because v[j] and v[k] are orthogonal, and the others are part of the definition of A[j](n).
2 (π^(2n+1))/(2n+1)
- 2A[j](n)A[j](n)B[j]^2 + A[j](n)^2B[j]^2
- 2A[k](n)A[k](n)B[k]^2 + A[k](n)^2B[k]^2
In other words
B[n]^2 =2 (π^(2n+1))/(2n+1) - sum ( A[j](n)^2B[j]^2 )
From here, we can construct everything that we need.
Now, to solve the problem of sin(x), we'll need to calculate the appropriate inner products. We know that sin(x) is an odd function, which cuts away half of our worries. And because we are dealing with inner products, we can simply calculate in a convenient basis, then apply the transform above. So we look at , , .
These are evaluated using the technique "integration by parts". My highschool text taught this using u's and dv's and generally turning my head about, but a remainder approach works just as well.
d/dx[- x^5 cos(x)] = x^5 sin(x) -5 x^4 cos(x)
d/dx[5 x^4 sin(x)] = 5 x^4 cos(x) +20 x^3 sin(x)
d/dx[20 x^3 cos(x)] = -20 x^3 sin(x) +60 x^2 cos(x)
d/dx[-60 x^2 sin(x)] = -60 x^2 cos(x) -120 x sin(x)
d/dx[-120 x cos(x)] = -120 x sin(x) -120 cos(x)
d/dx[120 sin(x)] = 120 cos(x)
Int[ x^5 sin(x) ] = ( 5 x^4 - 60 x^2 + 120) sin(x) - (x^5 - 20 x^3 + 120 x ) cos(x)
Of course, the sin(x) terms vanish at [-π,π],cos(π) flips the sign, and the result is odd so we can simply take half the integral and double it: 2 ( π^5 - 20 π ^3 + 120 π ).
This step done, we are through with exact answers - I started dumping things into eXcel at this point to track the calculations. But it does all work out, producing the answer in the book. Hooray for that.
Now, the answer in the text, 0.987682x - 0.155171x^3 + .00564321x^5 compares favorably with the Taylor series over the entire range - in the taylor series, the x^5 term becomes dominant too soon, producing a significant error at π. But the Taylor series is more accurate in the middle of the range, which is what you would care about if you were approximating infinitesmal angles.
Finally, what about cosine? The process for approximating the polynomial is the same as before, using the even basis vectors instead of the odds.
Now, I was rather startled when I started working through cosine when I found that the coefficient for [1] was zero! But it's ok - what that really means is that the integral of cosine over this range is zero, which is right. The x^0 terms we expect are really coming from v[2] and v[4]. When we grind it out, we get a polynomial that is not the derivative of the sin(x) series we just found.
And that's right - again, the Taylor series has the best approximation in the middle, but huge errors at the endpoints. The derivative gives a slightly worse approximation in the middle, and loses some accuracy at the endpoints. The calculated polynomial isn't quite as good in the middle, but makes it up at the endpoints, where the error is roughly one quarter of the sin(x) derivative.
December 14, 2003 1:45 PM
| TrackBack