Jacobi method unsucessful #202
Description
Describe the bug
When computing the eigenvalues and eigenvectors using vsl.la.jacobi, with known solutions, the method gives false eigenvalues.
The solutions are as 0.5, 1.5, 2.5, 3.5 and so on.
Expected Behavior
Expected correct eigenvalues, and their corresponding Eigen-vectors.
Current Behavior
Solutions are given wrong values
Reproduction Steps
module main
import vsl.la
import vsl.vlas
import os
const xmin = -6.0
const xmax = 6.0
fn potential(x f64) f64 {
return x * x
}
fn flatten(m &la.Matrix[f64]) []f64 {
mut flat := []f64{}
for i in 0 .. m.m {
flat << m.get_row(i)
}
return flat
}
enum Method {
three_point
five_point
invalid
}
/*
* Create a matrix for the 1D Schrodinger equation
@param[in] n The number of points
@param[in] ve The potential energy function
@return corresponding matrix
*/
fn create_matrix(n int, ve fn (f64) f64, m Method) &la.Matrix[f64] {
/*
The equation:latex
\[
1/2* ( -d^2/dz^2 + v(z) ) psi(z) = E / (hbar omega) psi(z)
\]
*/
mut hamiltonian := la.Matrix.new[f64](n, n)
dx := (xmax - xmin) / f64(n - 1)
println('dx: ${dx}')
mut x := xmin
if m == .three_point {
c_1 := -1.0 / (dx * dx) / 2.0
c1 := -1.0 / (dx * dx) / 2.0
for i in 0 .. n {
if i == 0 {
// We set [0, 1] seperatlly
hamiltonian.set(i, i + 1, c1)
} else if i == n - 1 {
// We set [n-1, n-2] seperatlly
hamiltonian.set(i, i - 1, c_1)
} else {
// we set the off-set diagonal elements
hamiltonian.set(i, i - 1, c_1)
hamiltonian.set(i, i + 1, c1)
}
hamiltonian.set(i, i, (2.0 / (dx * dx) + ve(x)) / 2.0)
x += dx
}
} else if m == .five_point {
c_2 := 1.0 / (12 * dx * dx) / 2.0
c_1 := -16.0 / (12 * dx * dx) / 2.0
c1 := 1.0 / (12 * dx * dx) / 2.0
c2 := -16.0 / (12 * dx * dx) / 2.0
for i in 0 .. n {
if i == 0 {
hamiltonian.set(i, i + 1, c1)
hamiltonian.set(i, i + 2, c2)
} else if i == 1 {
hamiltonian.set(i, i - 1, c_1)
hamiltonian.set(i, i + 1, c1)
hamiltonian.set(i, i + 2, c2)
} else if i == n - 2 {
hamiltonian.set(i, i - 2, c_2)
hamiltonian.set(i, i - 1, c_1)
hamiltonian.set(i, i + 1, c1)
} else if i == n - 1 {
hamiltonian.set(i, i - 2, c_2)
hamiltonian.set(i, i - 1, c_1)
} else {
hamiltonian.set(i, i - 2, c_2)
hamiltonian.set(i, i - 1, c_1)
hamiltonian.set(i, i + 1, c1)
hamiltonian.set(i, i + 2, c2)
}
hamiltonian.set(i, i, (30 / (12 * dx * dx) + ve(x)) / 2.0)
x += dx
}
} else {
panic('Invalid method')
}
return hamiltonian
}
fn main() {
n := 20
mut matrix := create_matrix(n, potential, Method.three_point)
mut eigen_vectors := la.Matrix.new[f64](n, n)
mut values := []f64{len: n}
la.jacobi(mut eigen_vectors, mut values, mut matrix) or { panic(err) } // TODO: This method does not work?
println(values)
Possible Solution
Unknown
Additional Information/Context
The matrix
The matrix in question is composed of a 3 or five-point stensil with a quadratic contribution to the main diagonal.
Thus, the main diagonal elements change, but the matrix is symmetric and banded, with either bandwidth 3 or five depending on stencil, both implemented in the code snippet above
Other efforts
Using vsl.vlas.dgeev was also unsuccessful, possibly due to lack of documentation, also in vsl.vlas.lapack_common.v C.LAPACKE_dsyev has implementation native to the package so that solution method is also rendered not possible.
V version
V 0.4.5
Version used
Latest ( v upgrade ran before submitted)
Environment details (OS name and version, etc.)
V full version: V 0.4.5 4bc9a8f.d692d88
OS: macos, macOS, 14.0, 23A344
Processor: 8 cpus, 64bit, little endian, Apple M2
getwd: /Users/andreas/Desktop/School/fk8029/report3/code
vexe: /Users/andreas/v/v
vexe mtime: 2024-04-07 13:53:21
vroot: OK, value: /Users/andreas/v
VMODULES: OK, value: /Users/andreas/.vmodules
VTMP: OK, value: /tmp/v_501
Git version: git version 2.39.3 (Apple Git-145)
Git vroot status: weekly.2024.09-208-gd692d884 (8 commit(s) behind V master)
.git/config present: true
CC version: Apple clang version 15.0.0 (clang-1500.0.40.1)
thirdparty/tcc status: thirdparty-macos-arm64 a668e5a0