Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jacobi method unsucessful #202

Open
thesombady opened this issue Apr 8, 2024 · 0 comments
Open

Jacobi method unsucessful #202

thesombady opened this issue Apr 8, 2024 · 0 comments
Labels
bug Something isn't working

Comments

@thesombady
Copy link

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

@thesombady thesombady added the bug Something isn't working label Apr 8, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

1 participant