How to define the Galerkin method for a general case?

[leocastel]

Thank you very much for providing such a nice Julia's package.
I would like to know how can I calculate the matrices in the Galerkin method for a general case.
For example, how can I calculate the following term:

$$L_{ij}=\int_a^b\frac{\phi_i(x)}{x^2}\frac{d\phi_j(x)}{dx} dx,$$

[jipolanco]

That's a good question, the computation of that kind of terms is not currently implemented, but it shouldn't be too difficult to add.

Currently the computation of matrices such as

$$M_{ij}=\int_a^b\phi_i(x) \frac{d\phi_j(x)}{dx} dx,$$

is done by adding the integrals computed in-between spline knots $t_i$. Roughly speaking:

$$M_{ij} = \sum_{i = 1}^N \int_{t_i}^{t_{i + 1}} \phi_i(x)\frac{d\phi_j(x)}{dx} dx$$

Since, in each $[t_i, t_{i + 1}]$ segment a B-spline of order $k$ is a polynomial of degree $n = k - 1$, then the integrand is (at most) a polynomial of degree $2n$ (well, actually $2n - 1$ in this case since we have one derivative). In that case, using an $n$-point Gauss-Legendre quadrature gives the exact integral (up to roundoff error).

In more general cases as in your example, Gauss-Legendre quadratures won't give the exact result (since the integrand is no longer a polynomial), but I guess they can still give you a good approximation.

You can take a look at the implementation. To modify the integrand, one would simply modify the following line:

@inbounds A[i, j] += y * bi * bj

by something like

@inbounds A[i, j] += y * bi * bj / x^2

I will try to add a more generic implementation (also allowing to change the size of the quadrature) when I find the time.

[leocastel]

Thank you for your fast reply.

So, the problem is related to the quadrature procedure that might have errors in the integration.
Perhaps, I can calculate the general terms by using another integration method. I will try this approach and let you know about it.

Best,

[jipolanco]

Update: this will be easier to do after #82 is merged.

In your example, you would do:

using BSplineKit
B = BSplineBasis(...)
f(x) = 1/x^2
L = galerkin_matrix(f, B, (Derivative(0), Derivative(1)))

One can also change the number of quadrature nodes to make sure that things converge. For example:

L = galerkin_matrix(f, B, (Derivative(0), Derivative(1)); quadrature = Galerkin.gausslegendre(Val(8)))

But l me know if you have any comments or if you find another approach.

[leocastel]

I am trying to solve the 2D Hydrogen atom using B-Splines, therefore I have to solve the following differential equation:

$$\left[\frac{d^2}{dr^2}+\frac{1}{r}\frac{d}{dr}+\frac{2}{r}-\frac{m^2}{r^2}\right]\psi_m(r)=-\epsilon \psi_m(r).$$

If I expand

$$\psi_m(r)=\sum_j c_i B_j(r)$$

in B-Splines, I would get the following generalized eigen problem:

$$Lc=-\epsilon R c,$$

where L and R are the Galerkin matrices for the lhs and rhs, respectively.

I tried to use your last suggestion, but I believe the result has some problems.
I've got a matrix for the lhs that the eigen command from LinearAlgebra package does not work.

I will another approach.
Thanks for the help.
Cheers,

[leocastel]

Interestingly, if I solve the simple problem:

$$\frac{d^2}{dr^2}\psi(r)=-\epsilon \psi(r),$$

I got the right answer for the following implementation of the second derivative:

Rgal = galerkin_matrix(R, (Derivative(1), Derivative(1))).

But it doesn't work if I use the following implementation of the second derivative:

Rgal = galerkin_matrix(R, (Derivative(2), Derivative(0))).

Here is the code:
``
knots_in = range(0, 1; length = 10)
using BSplineKit
using LinearAlgebra
B = BSplineBasis(BSplineOrder(4), knots_in)
using CairoMakie
using LaTeXStrings
CairoMakie.activate!(type = "svg")
R = RecombinedBSplineBasis(B, Derivative(0))
Sgal = galerkin_matrix(R)
Rgal = galerkin_matrix(R, (Derivative(1), Derivative(1)))
H=eigen(Rgal,Sgal)
vals, vecs = H
println(vals)
plot(vecs[:,1])

[jipolanco]

Can you check the types of the two Rgal versions? The first one should be a Hermitian matrix (because it's guaranteed to be symmetric), and maybe that's why it works...

In that case I'd guess it's a problem in BandedMatrices.jl, which only implements eigen(A, B) for A::Hermitian. Can you try converting the Galerkin matrix to a regular dense Matrix or to a SparseMatrixCSC from SparseArrays?

[leocastel]

You are completely correct. The first matrix is Hermitian, while the second is not.
In theory, both matrices should be Hermitian. I believe there is a problem with the integration method, right?

I don't know how to use the conversion SparseMatrixCSC you suggested.

In your implementation of the Gauss-Legendre method the range is always [-1,1] or do you rescale it for a different range of points?
Also, I would like to ask you how can I get the Bspline at a point x within the range.
Perhaps, I can try to investigate the integration and get Hermitian matrices.

[jipolanco]

I believe there is a problem with the integration method, right?

No, it's simply that galerkin_matrix tries to determine at compile time whether the matrix is Hermitian. So, if we can't know a priori if the matrix will be Hermitian or not, then the returned matrix will not be classified as Hermitian. This is for type stability reasons (we want to be sure that, given the types of the inputs, the output has one unique type). Right now, galerkin_matrix assumes that if both Derivative arguments are the same, then the matrix is of Hermitian type. Otherwise, the matrix may not be Hermitian, so it's treated as a regular banded matrix.

Note that here I'm only talking about the Hermitian type in Julia. One can have a matrix that is not Hermitian (the type) but which is Hermitian (the mathematical property). For example (adapted from the docs of ishermitian):

julia> a = [1 2; 2 -1]
2×2 Matrix{Int64}:
 1   2
 2  -1

julia> ishermitian(a)
true

julia> a isa Hermitian
false

If you know that the matrix is actually Hermitian, you can manually do:

Rgal = galerkin_matrix(R, (Derivative(2), Derivative(0))) 
@assert !(Rgal isa Hermitian)  # not "classified" as Hermitian, even though it actually is
@assert ishermitian(Rgal)      # check that mathematically, the matrix is actually Hermitian
Rgal_herm = Hermitian(Rgal)    # this should make `eigen` happy
H = eigen(Rgal_herm, Sgal)

I don't know how to use the conversion SparseMatrixCSC you suggested.

You should be able to do something like:

using SparseArrays
Rgal_sparse = sparse(Rgal)
Sgal_sparse = sparse(Sgal)
H = eigen(Rgal_sparse, Sgal_sparse)

I haven't checked whether eigen works in this case though.

In your implementation of the Gauss-Legendre method the range is always [-1,1] or do you rescale it for a different range of points?

It is translated and rescaled from $[-1, 1]$ to the range $[t_i, t_{i + 1}]$ mentioned in my original answer.

[leocastel]

Thank you very much for your reply.
Now, the simple problem:

$$\frac{d^2}{dr^2}\psi(r)=-\epsilon \psi(r),$$

works for both implementations of the second derivative:

Rgal = galerkin_matrix(R, (Derivative(1), Derivative(1))),

and

Rgal = galerkin_matrix(R, (Derivative(2), Derivative(0))).

Concerning the integration method, I was not completely clear when I wrote the following:

I believe there is a problem with the integration method, right?

I was referring to the case where there are non-polynomial terms such as,
f(x) = 1/x; R1= galerkin_matrix(f, R, (Derivative(1), Derivative(0)); quadrature = Galerkin.gausslegendre(Val(10))) g(x) = 2/x; R2= galerkin_matrix(g, R, (Derivative(0), Derivative(0)); quadrature = Galerkin.gausslegendre(Val(10)))
In these cases, there will be terms proportional to 1/x in the integral that could be wrongly calculated by the Gauss-Legendre quadrature.

The approach mentioned by you in the last reply,

using SparseArrays
Rgal_sparse = sparse(Rgal)
Sgal_sparse = sparse(Sgal)
H = eigen(Rgal_sparse, Sgal_sparse)

doesn't work because of the command eigen, was you have already mentioned.

[jipolanco]

I was referring to the case where there are non-polynomial terms such as,
f(x) = 1/x; R1= galerkin_matrix(f, R, (Derivative(1), Derivative(0)); quadrature = Galerkin.gausslegendre(Val(10))) g(x) = 2/x; R2= galerkin_matrix(g, R, (Derivative(0), Derivative(0)); quadrature = Galerkin.gausslegendre(Val(10)))
In these cases, there will be terms proportional to 1/x in the integral that could be wrongly calculated by the Gauss-Legendre quadrature.

Right, sorry. Yes, in that case Gauss-Legendre quadratures won't be exact any more. In general, I'd still expect them to provide a good approximation and to converge to the actual value as you increase the number of quadrature nodes. But this may not be the case if the singularity at $x = 0$ is in your integration interval.

[leocastel]

First of all, I'm sorry for my lack of knowledge about B-Splines at the beginning of my discussions.
Now, I have implemented a code from scratch that implements the B-Splines and the Garlekin matrices, which is in perfect agreement with your code.
I found out today that the term related to the Laplacian in polar coordinates that results in the matrix element of the Garlekin's matrix:

$$\int B_i(r)\frac{1}{r}\frac{d}{dr}B_j(r)\neq\int B_j(r)\frac{1}{r}\frac{d}{dr}B_i(r),$$

is the problematic one. The reason is that the Garlekin matrix will not be Hermitian for this term and therefore the eigenvalues are complex. Perhaps, this problem can be circumvented if I use the collocation point method. What do you think about it?
Best regards,

[leocastel]

Today, I implemented the collocation points method and finally got the expected results.
Thank you for your help.