Jacobi method unsucessful

[thesombady]

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