- #1
- 1,388
- 12
I have an integral equation of the form
y(x) = \int_0^{2\pi} dx'~K(x,x';\omega,\alpha)y(x'),
where K(x,x';\omega,\alpha) is a real, non-symmetric, non-separable kernel that depends on the parameters \omega and \alpha, as well as some others that I won't need to vary as much. The kernel is, unfortunately, relatively ugly, and I can't say too much about its properties analytically.
I have solved the equation numerically (in C) using the usual methods: I discretized the integral using Gauss-Legendre weights, w_j, on the [0,2*pi) interval, and then found the eigenvalues and eigenvectors of the matrix K_{ij} = K(x_i,x_j:\omega,\alpha)w_j, and picked the eigenvector corresponding to the eigenvalue of 1 (or close to it), which appears to be the largest eigenvector. I used the GSL routines to compute the weights and solve the eigensystem. With the corresponding eigenvector in hand, the interpolating formula for the solution is
y(x) = \sum_{j=0}^{N-1} K(x,x_j)w_j y_j,
for an N-point discretization, with y_j the eigenvector corresponding to an eigenvalue of 1.
However, there are a few problems with this approach:
(1) As \omega becomes small (less than ~100 for typical parameter choices I have been using), the kernel becomes increasingly sparse, being significantly non-zero only near (0,2Pi), (2Pi,0) and a curve near the diagonal. As a result, for small omega I need to sample very many points when I discretize the integral, otherwise the peaked regions are not well resolved and solving the eigensystem gives a largest eigenvalue that isn't even close to 1. For \omega greater than about 100, the kernel is not very sparse and ~30 discretization points is sufficient to get the right eigenvalue. Near a value of about 10, I need around 100 discretization points to get an eigenvalue near 1. Below a value of 10, the number of points I need is simply too large.
(2)The other problem is that I need to change the parameters \omega and \alpha quite frequently, which means that every time I change the parameters I need to regenerate the discretized kernel and resolve the eigensystem. A lot of the interesting features I am trying to look at also occur for values of omega below 100, which means that the matrices are typically of size ~100x100. The GSL routines I am using will compute all the eigenvalues and eigenvectors associated with the discrete kernel, but I only need the one eigenvector with eigenvalue (approximately) 1, which seems to be the largest and only real eigenvalue. There doesn't seem to be an option to have the routine just compute the largest eigenvalue and corresponding eigenvector (and I'm not sure if it would be significantly faster not to compute the others, or if it just uses less memory)
So, what I am looking for are methods to quickly compute the solution to an integral equation, which could be done many times over.
I would be interested in references to non-matrix methods that could solve the problem more quickly than the matrix methods, (and hopefully give me an interpolating formula for the result), or if I am to use matrix methods, c packages that can quickly compute the largest eigenvalue and eigenvector, so as to offset the need for many more points when \omega becomes small, or any other tricks that might speed up computation but won't take another month to implement.
Thanks for any help,
--Mute
y(x) = \int_0^{2\pi} dx'~K(x,x';\omega,\alpha)y(x'),
where K(x,x';\omega,\alpha) is a real, non-symmetric, non-separable kernel that depends on the parameters \omega and \alpha, as well as some others that I won't need to vary as much. The kernel is, unfortunately, relatively ugly, and I can't say too much about its properties analytically.
I have solved the equation numerically (in C) using the usual methods: I discretized the integral using Gauss-Legendre weights, w_j, on the [0,2*pi) interval, and then found the eigenvalues and eigenvectors of the matrix K_{ij} = K(x_i,x_j:\omega,\alpha)w_j, and picked the eigenvector corresponding to the eigenvalue of 1 (or close to it), which appears to be the largest eigenvector. I used the GSL routines to compute the weights and solve the eigensystem. With the corresponding eigenvector in hand, the interpolating formula for the solution is
y(x) = \sum_{j=0}^{N-1} K(x,x_j)w_j y_j,
for an N-point discretization, with y_j the eigenvector corresponding to an eigenvalue of 1.
However, there are a few problems with this approach:
(1) As \omega becomes small (less than ~100 for typical parameter choices I have been using), the kernel becomes increasingly sparse, being significantly non-zero only near (0,2Pi), (2Pi,0) and a curve near the diagonal. As a result, for small omega I need to sample very many points when I discretize the integral, otherwise the peaked regions are not well resolved and solving the eigensystem gives a largest eigenvalue that isn't even close to 1. For \omega greater than about 100, the kernel is not very sparse and ~30 discretization points is sufficient to get the right eigenvalue. Near a value of about 10, I need around 100 discretization points to get an eigenvalue near 1. Below a value of 10, the number of points I need is simply too large.
(2)The other problem is that I need to change the parameters \omega and \alpha quite frequently, which means that every time I change the parameters I need to regenerate the discretized kernel and resolve the eigensystem. A lot of the interesting features I am trying to look at also occur for values of omega below 100, which means that the matrices are typically of size ~100x100. The GSL routines I am using will compute all the eigenvalues and eigenvectors associated with the discrete kernel, but I only need the one eigenvector with eigenvalue (approximately) 1, which seems to be the largest and only real eigenvalue. There doesn't seem to be an option to have the routine just compute the largest eigenvalue and corresponding eigenvector (and I'm not sure if it would be significantly faster not to compute the others, or if it just uses less memory)
So, what I am looking for are methods to quickly compute the solution to an integral equation, which could be done many times over.
I would be interested in references to non-matrix methods that could solve the problem more quickly than the matrix methods, (and hopefully give me an interpolating formula for the result), or if I am to use matrix methods, c packages that can quickly compute the largest eigenvalue and eigenvector, so as to offset the need for many more points when \omega becomes small, or any other tricks that might speed up computation but won't take another month to implement.
Thanks for any help,
--Mute
