I'm trying to reproduce the sin(x) approximation example from "Linear algebra done right" using Maxima, but I get weird results. Here is my code:
load("eigen");
ip (f, g) := integrate (f * g, u, a, b);
basis: makelist (u^k, k, 0, 5);
e: gramschmidt (basis, ip), a=-%pi, b=%pi;
sl: makelist (sin(u), x, 1, 6);
re: map(ip, sl, e)*e, a=-%pi, b=%pi;
expand(float(apply("+", re)));
which gives me 4.867477078388566*u^5 - 66.27194960446585*u^3 + 195.5285758118965*u
while the correct result is 0.987862*x − 0.155271*x^3 + 0.00564312*x^5
What am I doing wrong?
Looks like you need to normalize the basis before applying it to sin(u). Here's how I phrased it; no doubt there are other ways to say it.
(%i2) load("eigen") $
(%i3) ip (f, g) := integrate (f * g, u, -%pi, %pi) $
(%i4) n: 5 $
(%i5) basis: makelist (u^k, k, 0, n);
2 3 4 5
(%o5) [1, u, u , u , u , u ]
(%i6) e: gramschmidt (basis, ip);
2 2 2 2
3 u - %pi u (5 u - 3 %pi )
(%o6) [1, u, ───────────, ─────────────────,
3 5
4 2 2 4 4 2 2 4
35 u - 30 %pi u + 3 %pi u (63 u - 70 %pi u + 15 %pi )
───────────────────────────, ────────────────────────────────]
35 63
Find normalizing factors z and use them to construct an orthonormal basis.
(%i7) z: map (lambda ([ee], ip (ee, ee)), e);
3 5 7 9 11
2 %pi 8 %pi 8 %pi 128 %pi 128 %pi
(%o7) [2 %pi, ──────, ──────, ──────, ────────, ─────────]
3 45 175 11025 43659
(%i8) e1: e / sqrt(z);
2 2
1 sqrt(3) u sqrt(5) (3 u - %pi )
(%o8) [─────────────────, ──────────────, ─────────────────────,
sqrt(2) sqrt(%pi) 3/2 3/2 5/2
sqrt(2) %pi 2 %pi
2 2 4 2 2 4
sqrt(7) u (5 u - 3 %pi ) 3 (35 u - 30 %pi u + 3 %pi )
─────────────────────────, ───────────────────────────────,
3/2 7/2 7/2 9/2
2 %pi 2 %pi
4 2 2 4
sqrt(11) u (63 u - 70 %pi u + 15 %pi )
─────────────────────────────────────────]
7/2 11/2
2 %pi
Verify that e1 is an orthonormal basis.
(%i9) genmatrix (lambda ([i, j], ip (e1[i], e1[j])), n + 1, n + 1);
┌ ┐
│ 1 0 0 0 0 0 │
│ │
│ 0 1 0 0 0 0 │
│ │
│ 0 0 1 0 0 0 │
(%o9) │ │
│ 0 0 0 1 0 0 │
│ │
│ 0 0 0 0 1 0 │
│ │
│ 0 0 0 0 0 1 │
└ ┘
Now use the basis to construct an approximation to sin(u).
(%i10) coeffs: map (lambda ([ee], ip (sin(u), ee)), e1);
3
sqrt(2) sqrt(3) sqrt(7) (4 %pi - 60 %pi)
(%o10) [0, ───────────────, 0, ─────────────────────────, 0,
sqrt(%pi) 3/2 7/2
2 %pi
5 3
sqrt(11) (16 %pi - 1680 %pi + 15120 %pi)
──────────────────────────────────────────]
7/2 11/2
2 %pi
(%i11) approx: coeffs . e1;
5 3
(%o11) (11 (16 %pi - 1680 %pi + 15120 %pi) u
4 2 2 4 11
(63 u - 70 %pi u + 15 %pi ))/(128 %pi )
3 2 2
7 (4 %pi - 60 %pi) u (5 u - 3 %pi ) 3 u
+ ───────────────────────────────────── + ────
7 2
8 %pi %pi
(%i12) expand (approx);
5 5 5 3 3
693 u 72765 u 654885 u 315 u 39375 u
(%o12) ────── - ──────── + ───────── - ────── + ────────
6 8 10 4 6
8 %pi 8 %pi 8 %pi 4 %pi 4 %pi
3
363825 u 105 u 16065 u 155925 u
- ───────── + ────── - ─────── + ────────
8 2 4 6
4 %pi 8 %pi 8 %pi 8 %pi
(%i13) fpprintprec: 6;
(%o13) 6
(%i14) float (expand (approx));
5 3
(%o14) 0.00564312 u - 0.155271 u + 0.987862 u