Orthonormalizing Functions on a Given Interval: A MATLAB Approach

[member 428835]

Hi PF!

I recently wrote a program in MATLAB to orthonormalize continuous functions on a given interval $x\in[a,b]$. Let's call the initial linearly independent functions I seek to orthonormalize $L_i$ and let's call these functions $\psi_i$ after I orthonormalize them.

After running my script, I perform an integration check, where I create an $N\times N$ matrix $A:A_{ij} = \int_a^b \psi_i\psi_j\,dx$, $N$ being the number of $L_i$'s I use. Theoretically $A$ should be the identity matrix.

When I let $N=5$, my maximum error outside the diagonals for a particular sequence of functions (too difficult to explain which functions) is very small $O(10^{-5})$. However, when I let $N=20$, the maximum error off the diagonals is $O(1)$.

My question is, do any of you know a set of functions that, once orthonormalized, give another analytic set of functions? I thought about taking $L_i = x^i:x\in[-1,1]$ to generate the Legendre polynomials, but these polynomials are not orthonormal on $x\in[-1,1]$. I used $\sin(nx)$ and $\cos(nx)$ and they work, but not all sequences of functions seem to (like the first function I eluded to in the previous paragraph).

Does anyone know either what the Legendre polynomials are orthonormal over OR a set of linearly independent functions that have a known analytic corresponding orthonormal set of functions?

I'm trying to check my code and this would be very helpful. For your info, here's the script I wrote

N = 20;% number of functions

dx = 0.001;
x = -pi:dx:pi;
Li = zeros(N,size(x,2));
for n = 1:size(Li,1)
    Li(n,:) = cos(n*x)./sqrt(pi);% linearly independent functions to orthonormalize
end

psi = zeros(size(Li));
ortho = zeros(size(psi));

% gram-schmidt
for ii = 1:size(psi,1)
        for s = 1:ii-1
            ortho(ii,:) = ortho(ii,:)+psi(s,:)*trapz(psi(s,:).*Li(ii,:))*dx;
        end% for s
        psi(ii,:) = (Li(ii,:)-ortho(ii,:))./sqrt(trapz((Li(ii,:)-ortho(ii,:)).^2)*dx);% Gram-Schmidt on Li to construct psi
end% for ii

 % check orthogonality

for kk = 1:size(Li,1)
    for jj = 1:size(Li,1)
        orBC(jj,kk) = trapz(psi(kk,:).*psi(jj,:))*dx;% orthogonality check on psi (should be Identity matrix)
    end% for jj
end% for kk

for kk = 1:size(Li,1)
    for jj = 1:size(Li,1)
        orBCexact(jj,kk) = trapz(Li(kk,:).*Li(jj,:))*dx;% orthogonality check on Li (should be Identity matrix)
    end% for jj
end% for kk

 

[NFuller]

Does anyone know either what the Legendre polynomials are orthonormal over OR a set of linearly independent functions that have a known analytic corresponding orthonormal set of functions?

I'm not completely sure what you are asking but the Legendre polynomials are orthogonal over the interval $[-1,1]$ such that
$$\int_{-1}^{1}P_{n}(x)P_{m}(x)dx=\frac{2}{2n+1}\delta_{nm}$$

 

[member 428835]

I'm not completely sure what you are asking but the Legendre polynomials are orthogonal over the interval $[-1,1]$ such that
$$\int_{-1}^{1}P_{n}(x)P_{m}(x)dx=\frac{2}{2n+1}\delta_{nm}$$

Thanks for replying! Yea, I saw that orthonormality relationship. So after closer inspection I think the script I attached works for most functions (it certainly does for Legendre polynomials. However, I have an interesting set of functions, which are linear combinations of Legendre polynomials over $[-1,1]$, which, when I enter them in as $L_i$, they are not all orthogonal.

Let me explain my problem in detail and if you don't mind you can help me out. Let's call the $ith$ Legendre polynomial $p_i:i=0,1,2...$. I have a piecewise function $f$ defined as
$$
f = \begin{cases}
f_1 = \sum_{i=0}^N b_i p_i(x)&x\in [-1,\zeta_1] \\
f_3 = 0 & x \in (\zeta_1,\zeta_2)\\
f_2 = \sum_{i=0}^N c_i p_i(x)&x\in [\zeta_2,1]
\end{cases}
$$
I have three boundary conditions $f$ is subject to. In matrix form, these three boundary conditions form a $3\times2(N+1)$ matrix $M$ that operates on the vector coefficients $\mathbf{b} = [b_0,b_1...b_N]^T$ and similarly for $\mathbf{c}$. In matrix form this looks like $$
[M]\begin{bmatrix}
\mathbf{b}\\
\mathbf{c}
\end{bmatrix} = \mathbf{0}.
$$
To solve this I compute $null(M)$, which is a $2(N+1)\times 2N-1$ matrix. Then there are $2N-1$ coefficient vectors $[\mathbf{b}; \mathbf{c}]^T$ that satisfy the boundary condition matrix $M$. So, once computing $null(M)$, I know $\mathbf{b}, \mathbf{c}$, and thus I know $f_1$ and $f_2$ and consequently I know $f$. Since there are $2N-1$ coefficient vectors, I actually have $2N-1$ linearly independent $f$ functions. However, these functions $f$ are not necessarily orthonormal over $[-1,1]$. This is where I apply Gram-Schmidt. When I apply Gram-Schmidt, it will transform $f_i\to\psi_i:i=0,1,2...2N-1$ where $\psi_i$ should be orthonormal.

To verify each $\psi_i$ is othonormal over $x\in[-1,1]$, I make a matrix $A : A_{ij} = \int_{-1}^1 \psi_i \psi_j\,dx$. If everything is correct, $A = I_{2N-1}$, the identity matrix. When $N = 5$ everything works out very well, and $A \approx I$. However, when $N=10$, $A \neq I$, where early entries, say $A_{3,4}$ or $A_{2,7}$ are approximately zero yet higher entries $A_{17,18}\approx 1$. I should say $A_{ii}=1$ for any $N$ selected.

When troubleshooting this, I was curious about the vector coefficients from $null(M)$. So I plotted $\psi_{15-19}$ and noticed they are almost identical. Any idea what's going on here? Attached are three plots: $\psi_i:i=0,1,...,19$, $\psi_i:i=15,...,19$, and a zoomed in view of $\psi_i:i=15,...,19$, all plotted over $x$.

Sorry, I keep forgetting how to post a picture directly on this post.

 

[jasonRF]

EDIT: I was looking at your code again, and it looks like may be already using the modified Gram-Schmidt? If so then my first guess is wrong. But to be honest I am confusing myself trying to look at your code - a problem with my brain at the end of a long day, not your code!

One issue that might be causing you grief is that, even for the "simple" case of finite vectors (not functions) the classical Gram-Schmidt algorithm is numerically unstable on a floating point machine. You are essentially implementing classical Gram-Schmidt so that could be the problem; your use of trapz for estimation of the inner products may also be adding to the instability? Anyway, see the Wikipedia page:

https://en.wikipedia.org/wiki/Gram–Schmidt_process

Your symptom of having higher order $\psi$ functions being very similar is the same thing I have seen before when Gram-Schmidt is unstable. You might would not see the instability if you start with a nice set of functions on a nice interval (say, sinusoids on a [-pi,pi] or Legendre polynomials on [-1,1]) except for very large N, but for some choices of functions you will see it for N relatively small, like 10.

jason

 

[member 428835]

This is exactly the advice I was looking for. I'll do some research, thanks!

 

[member 428835]

So I used the algorithm attached on wikipedia, and for most functions I get the same result as my code in post 1. However, when I apply this to the function I refer to in post 1, (which I call $f$ in the code below) the plot looks nasty for large $\psi$ functions. The largest value of $A$ off the diagonal is 0.19, so that is better, but still not small enough to consider it orthogonal. Attached is a picture of the new $\psi$ values for the function.

Any idea what's happening? My code in post 1 seems to be more stable, as the "orthonormal" functions appear to be continuous, yet obviously my code gives linearly dependent solutions for large $N$. What's going on here? I'm happy to provide any additional info that would help.

I could program this in Mathematica and avoid all numerical errors, since the function $f$ that's giving me trouble is a linear combination of Bessell functions, but it would be so nice to get this working in MATLAB. Plus, this is a good learning experience for me.

V = f';
  n = size(V,1);
    k = size(V,2);
    U = zeros(n,k);
    U(:,1) = V(:,1)/sqrt(trapz(V(:,1).*V(:,1))*dx);
    for i = 2:k
      U(:,i) = V(:,i);
      for j = 1:i-1
        U(:,i) = U(:,i) - ( U(:,i)'*U(:,j) )/( U(:,j)'*U(:,j) )*U(:,j);
      end
      U(:,i) = U(:,i)/sqrt(trapz(U(:,i).*U(:,i))*dx);
    end
U = U';

 

[jasonRF]

One thing about that code is that you are not consistent about when you use a trapz inner product (defining U(:,1) and normalizing U(:,i) at the end of each for i= ... loop) versus a 'normal' inner product (inside the for j= ... loop).

Have you tried running this using simple inner products first instead of your trapz inner products?

jason

 

[member 428835]

One thing about that code is that you are not consistent about when you use a trapz inner product (defining U(:,1) and normalizing U(:,i) at the end of each for i= ... loop) versus a 'normal' inner product (inside the for j= ... loop).

Have you tried running this using simple inner products first instead of your trapz inner products?

jason

Shoot, good call. So I've fixed the Gram-Schmidt as shown below but it's still giving me very discontinuous outputs nearly identical to those shown in post 6. I did run this on simpler Legendre polynomials and it works out well. Really confused why it's not working out here. I could send you a plot of the functions that are not orthonormalizing under this procedure if you'd like?

V = f';
  n = size(V,1);
    k = size(V,2);
    U = zeros(n,k);
    U(:,1) = V(:,1)/sqrt(trapz(V(:,1).*V(:,1))*dx);
    for i = 2:k
      U(:,i) = V(:,i);
      for j = 1:i-1   
         U(:,i) = U(:,i) - (trapz(U(:,i).*U(:,j))*dx )/( trapz(U(:,j).*U(:,j))*dx )*U(:,j);
      end
      U(:,i) = U(:,i)/sqrt(trapz(U(:,i).*U(:,i))*dx);
    end
U = U';

 

[jasonRF]

I'm looking at your code and that is definitely not the modified Gram-Schmidt - it is the unstable classical Gram-Schmidt. My fault for not making sure wikipedia was correct!

The stable modified Gram-Schmidt (here coded up as a QR decomposition) would look more like,

% QR decomposition via modified Gram-Schmidt.  V = Q*R
[m,n] = size(V);  % columns of V have initial vectors
R = zeros(n,n);   % R is upper triangular matrix
Q = zeros(m,n);   % columns of Q will have orthonormal vectors
for k=1:n
   R(k,k) = norm(V(:,k));
   Q(:,k) = V(:,k)/R(k,k);
   for jj=k+1:n
      R(k,jj) = Q(:,k)'*V(:,jj);
     V(:,jj) = V(:,jj) - Q(:,k)*R(k,jj);
   end
end

A trapz version would only have to modify two lines.

In any case, I recommend looking at your vectors with standard tools and standard (not trapz) inner products to make sure your problem is well conditioned. One way is to use Matlab's built-in QR decomposition. I suspect it does not use Gram-Schmidt in any form, but instead uses either Householder reflectors or Givens rotations, which are even better than the modified G-S, especially if orthogonality of your vectors is of prime concern. For your problem try,

% let columns of V be vectors to make orthonormal
[Q,R] = qr(V,0);  % columns of Q give orthonormalized set of vectors

If you try this and the matrix Q is orthogonal (with respect to standard inner products) to an acceptable degree, then try to use the modified Gram-Schmidt above and see how it does. If it is fine, then a trapz version of the modified Gram-Schmidt code should be adequate for you. If not, then you will need to decide if the trapz stuff is really that much better than a standard inner product. If it is, then you will need to code up something like a Householder algorithm with trapz, which will require you to re-derive parts because all of the versions you are likely to find will implicitly be assuming standard inner products. This is do-able, but will take some real work on your part.

Hope that helps.

Jason

EDIT: any book on numerical linear algebra should cover all of this. Examples are Golub and Van Loan, "Matrix Computations", and Watkins "Fundamentals of Matrix Computations". Any edition for either book.