lars/FYS4150
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Finish project 2

Browse Source
This commit is contained in:
Lars Bogner
2025-09-18 09:06:36 +02:00
parent e62c42bb2c
commit fededb4ed3
13 changed files with 688 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
.gitignore
+9 -1
View File
@@ -272,4 +272,12 @@ build/
src/project1/poisson
src/project1/poisson_solver
src/project1/poisson_timing
src/project1/rel_errors.txt
src/project1/rel_errors.txt
# Project 2 executables
src/project2/arma_eigenproblem_checker
src/project2/buckling_beam_solver
src/project2/off_diagonal_finder
src/project2/scaling_tester
src/project2/*.txt
projects/project2/main.pdf
BIN
View File
Binary file not shown.
projects/project2/main.tex
+182
View File
@@ -0,0 +1,182 @@
\documentclass[english,notitlepage]{revtex4-1} % defines the basic parameters of the document
%For preview: skriv i terminal: latexmk -pdf -pvc filnavn
\let\tablenum\relax % to avoid clash with siunitx
% if you want a single-column, remove reprint
% allows special characters (including æøå)
\usepackage[utf8]{inputenc}
%\usepackage[english]{babel}
%% note that you may need to download some of these packages manually, it depends on your setup.
%% I recommend downloading TeXMaker, because it includes a large library of the most common packages.
\usepackage{amssymb} % mathematical symbols (physics imports amsmath)
\usepackage{amsmath}
\usepackage{graphicx} % include graphics such as plots
\usepackage{xcolor} % set colors
\usepackage{hyperref} % automagic cross-referencing (this is GODLIKE)
\usepackage{listings} % display code
\usepackage{subfigure} % imports a lot of cool and useful figure commands
\usepackage{float}
%\usepackage[section]{placeins}
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{subfigure}
\usepackage{tikz}
\usepackage{siunitx}
\usetikzlibrary{quantikz}
% defines the color of hyperref objects
% Blending two colors: blue!80!black = 80% blue and 20% black
\hypersetup{ % this is just my personal choice, feel free to change things
colorlinks,
linkcolor={red!50!black},
citecolor={blue!50!black},
urlcolor={blue!80!black}}
%% Defines the style of the programming listing
%% This is actually my personal template, go ahead and change stuff if you want
%% USEFUL LINKS:
%%
%% UiO LaTeX guides: https://www.mn.uio.no/ifi/tjenester/it/hjelp/latex/
%% mathematics: https://en.wikibooks.org/wiki/LaTeX/Mathematics
%% PHYSICS ! https://mirror.hmc.edu/ctan/macros/latex/contrib/physics/physics.pdf
%% the basics of Tikz: https://en.wikibooks.org/wiki/LaTeX/PGF/Tikz
%% all the colors!: https://en.wikibooks.org/wiki/LaTeX/Colors
%% how to draw tables: https://en.wikibooks.org/wiki/LaTeX/Tables
%% code listing styles: https://en.wikibooks.org/wiki/LaTeX/Source_Code_Listings
%% \includegraphics https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
%% learn more about figures https://en.wikibooks.org/wiki/LaTeX/Floats,_Figures_and_Captions
%% automagic bibliography: https://en.wikibooks.org/wiki/LaTeX/Bibliography_Management (this one is kinda difficult the first time)
%% REVTeX Guide: http://www.physics.csbsju.edu/370/papers/Journal_Style_Manuals/auguide4-1.pdf
%%
%% (this document is of class "revtex4-1", the REVTeX Guide explains how the class works)
%% CREATING THE .pdf FILE USING LINUX IN THE TERMINAL
%%
%% [terminal]$ pdflatex template.tex
%%
%% Run the command twice, always.
%% If you want to use \footnote, you need to run these commands (IN THIS SPECIFIC ORDER)
%%
%% [terminal]$ pdflatex template.tex
%% [terminal]$ bibtex template
%% [terminal]$ pdflatex template.tex
%% [terminal]$ pdflatex template.tex
%%
%% Don't ask me why, I don't know.
\newcommand{\tderiv}[2]{\frac{\mathrm{d} #1}{\mathrm{d} #2}} % defines a command for total derivatives
\newcommand{\pderiv}[2]{\frac{\partial #1}{\partial #2}} % defines a command for partial derivatives
\newcommand{\ptderiv}[3]{\frac{\mathrm{d}^{#3} #1}{\mathrm{d} #2^{#3}}} % defines a command for total derivatives of arbitrary power
\newcommand{\ppderiv}[3]{\frac{\partial^{#3} #1}{\partial #2^{#3}}} % defines a command for partial derivatives of arbitrary power
\begin{document}
\title{Solutions to Project 2} % self-explanatory
\author{Lars Bogner} % self-explanatory
\date{\today} % self-explanatory
\noaffiliation % ignore this, but keep it.
\maketitle
\textit{Code hosted at \url{https://github.uio.no/larsbog/FYS4150}}
\section*{Problem 1}
Using the rescaling of the positions $\hat x \equiv \frac{x}{L}$ we define the following relations:
\begin{equation}
u_x \equiv u, \quad u_{\hat x} \equiv \hat u, \quad \text{with} \quad u(x) = \hat u(\hat x) = u (L \hat x).
\end{equation}
Now we can rewrite the equation as
\begin{equation}
\begin{aligned}
\gamma \ptderiv{u(x)}{x}{2} &= - F u(x)\\
\gamma \ptderiv{u(L \hat x)}{(L \hat x)}{2} &= - F u(L \hat x).
\end{aligned}
\end{equation}
In this step we simply replaced $x$ by the equivalent expression $L \hat x$. Now we can use the chain rule to rewrite the left hand side:
\begin{equation}
\gamma \ppderiv{u(L \hat x)}{\hat x}{2} \underbrace{\ppderiv{\hat x}{(L \hat x)}{2}}_{=L^{-2}} = - F u(L \hat x).
\end{equation}
Rearranging the factors gives the final result:
\begin{equation}
\ptderiv{\hat u (\hat x)}{\hat x}{2} = - \frac{F L^2}{\gamma} \hat u(\hat x) \equiv - \lambda \hat u(\hat x).
\end{equation}
\section*{Problem 2}
The code is implemented in \texttt{src/project2/arma\_eigenproblem\_checker.cpp}. The reusable functions for the tridiagonal matrix, the analytical solution and a wrapper around the Armadillo eigenvalue solver are implemented in \texttt{./src/helpers.cpp}.
\section*{Problem 3}
\subsection*{Subproblem (a)}
The function with the signature \texttt{double max\_offdiag\_symmetric(const arma::mat\& A, int\& k, int\& l)} finds the largest off-diagonal element in the symmetric matrix \texttt{A} and returns its value. The algorithm works like described in Algorithm \ref{alg:maxoffdiag}. The code is implemented in \texttt{./src/helpers.cpp}.
\begin{figure}
\begin{algorithm}[H]
\caption{Finding the largest off-diagonal element in a symmetric matrix}
\label{alg:maxoffdiag}
\begin{algorithmic}
\Function{max\_offdiag\_symmetric}{$A, k, l$}
\State $n \gets$ number of rows of $A$
\State $\mathrm{max} \gets 0$
\For{$i = 0$ to $n-1$}
\For{$j = i+1$ to $n-1$}
\If{$|A_{ij}| > \mathrm{max}$}
\State $\mathrm{max} \gets |A_{ij}|$
\State $k \gets i$
\State $l \gets j$
\EndIf
\EndFor
\EndFor
\State \Return $\mathrm{max}$
\EndFunction
\end{algorithmic}
\end{algorithm}
\end{figure}
\subsection*{Subproblem (b)}
The code correctly identifies the largest off-diagonal element as $0.7$ at position $(1, 2)$. The function call is tested in \texttt{src/project2/off\_diagonal\_finder.cpp}.
\section*{Problem 4}
The Jacobi rotation algorithm is implemented as described in the lecture notes. The implementation is found in \texttt{src/helpers.cpp} as the function \texttt{void jacobi\_eigensolver(const arma::mat\& A, double eps, arma::vec\& eigenvalues, arma::mat\& eigenvectors, const int maxiter, int\& iterations, bool\& converged)}. It uses the \texttt{jacobi\_rotate} function to perform the actual rotation. The algorithm is tested in \texttt{src/\-project2/\-arma\_\-eigenproblem\_\-checker.cpp}. It modifies the eigenvalues and eigenvectors passed by reference. The function also returns the number of iterations used and whether the algorithm converged within the maximum number of iterations by modifying the variables passed by reference.
\section*{Problem 5}
\subsection*{Subproblem (a)}
The scaling of the number of iterations is tested in \texttt{src/project2/scaling\_tester.cpp}. The results are shown in \autoref{fig:scaling}. On the logarithmic scale, we see a trend that is roughly of the order $\sqrt{\log N}$.
\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{include/scaling_data_iterations.pdf}
\caption{Scaling of the number of iterations with the matrix size $N$. }
\label{fig:scaling}
\end{figure}
\subsection*{Subproblem (b)}
The scaling is also tested for a dense matrix, i.e. a matrix without the tridiagonal structure. The results show nearly no dependance on the density of the matrix. This is to be expected, since the Jacobi algorithm creates new non-zero elements in each iteration, quickly filling up the matrix before conerging to the solution.
\section*{Problem 6}
\begin{figure}
\centering
\includegraphics[width=\textwidth]{include/eigen_data_10_plot.pdf}
\caption{Comparison of the numerical and analytical solution of the buckling beam problem for $N=10$.}
\label{fig:buckling_beam_10}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\textwidth]{include/eigen_data_100_plot.pdf}
\caption{Comparison of the numerical and analytical solution of the buckling beam problem for $N=100$.}
\label{fig:buckling_beam_100}
\end{figure}
The numerical and analytical solution of the buckling beam problem is compared in \autoref{fig:buckling_beam_10} for $N=10$ and in \autoref{fig:buckling_beam_100} for $N=100$. The numerical solution is in very good agreement with the analytical solution, even for the small value of $N=10$. For $N=100$, the solutions are visually indistinguishable. Nevertheless, the absolute error between the solutions is displayed on the right-hand side of the figures. The code to generate the plots is found in \texttt{src/project2/python/buckling\_beam\_plotter.py}.
\bibliography{include/citations}
\end{document}
projects/project2/mainNotes.bib
+2
View File
@@ -0,0 +1,2 @@
@CONTROL{REVTEX41Control}
@CONTROL{apsrev41Control,author="08",editor="1",pages="0",title="",year="1"}
src/project2/arma_eigenproblem_checker.cpp
+60
View File
@@ -0,0 +1,60 @@
#include "include/helpers.hpp"
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
int main(int argc, char** argv) {
int N = 6;
if (argc > 1){ // Allow user to specify matrix size during call
N = std::stoi(argv[1]);
}
// Setup of the matrix
double h = 1.0 / (N + 1);
mat A = create_tridiagonal_matrix(N, h);
cout << "Tridiagonal matrix A:" << endl << A << endl;
// Using Armadillo's built-in solver for comparison
vec eigenvalues;
mat eigenvectors;
// Calling arma
eig_sym(eigenvalues, eigenvectors, A);
// Sorting and normalizing the results
normalize_columns(eigenvectors);
sort_eigenpairs(eigenvalues, eigenvectors);
cout << "Eigenvalues:" << endl << eigenvalues << endl;
cout << "Normalized eigenvectors:" << endl << eigenvectors << endl;
// Generate analytical solution
vec analytical_eigenvalues;
mat analytical_eigenvectors;
analytical_solution(N, h, analytical_eigenvalues, analytical_eigenvectors);
// Again: Sort the pairs of eigenvalues and eigenvectors by eigenvalue magnitude
normalize_columns(analytical_eigenvectors);
sort_eigenpairs(analytical_eigenvalues, analytical_eigenvectors);
cout << "Analytical Eigenvalues:" << endl << analytical_eigenvalues << endl;
cout << "Analytical Normalized Eigenvectors:" << endl << analytical_eigenvectors << endl;
// Using self-made Jacobi solver
// Setup variables to capture output
vec jacobi_eigenvalues;
mat jacobi_eigenvectors;
int iterations;
bool converged;
double eps = 1e-8;
int maxiter = 1e5;
// Calling the solver
jacobi_eigensolver(A, eps, jacobi_eigenvalues, jacobi_eigenvectors, maxiter, iterations, converged);
if (converged) {
cout << "Jacobi method converged in " << iterations << " iterations." << endl;
} else {
cout << "Jacobi method did not converge within the maximum number of iterations." << endl;
}
normalize_columns(jacobi_eigenvectors); // Normalize the eigenvectors
sort_eigenpairs(jacobi_eigenvalues, jacobi_eigenvectors);
cout << "Jacobi Eigenvalues:" << endl << jacobi_eigenvalues << endl;
cout << "Jacobi Normalized Eigenvectors:" << endl << jacobi_eigenvectors << endl;
return 0;
}
src/project2/buckling_beam_solver.cpp
+85
View File
@@ -0,0 +1,85 @@
#include <armadillo>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "include/helpers.hpp"
using namespace std;
using namespace arma;
void write_boundary_conditions(ofstream& outfile, int num_eigenvectors, int decimal_precision){
for(int i = 0; i < num_eigenvectors; i++){
outfile << setprecision(decimal_precision) << scientific << 0.0 << "," << 0.0; // Boundary condition at x=1
if (i < num_eigenvectors - 1) {
outfile << ",";
}
else {
outfile << endl;
}
}
}
int main(int argc, char** argv) {
int n = 10;
const int NUM_EIGENVECTORS = 3; // Number of eigenvectors to save
const int DECIMAL_PRECISION = 12; // Decimal precision for output
if (argc > 1){
n = std::stoi(argv[1]);
}
int N = n - 1; // Number of interior points
double h = 1.0 / (n); // Step size
mat A = create_tridiagonal_matrix(N, h);
// Generate analytical solution
vec analytical_eigenvalues;
mat analytical_eigenvectors;
analytical_solution(N, h, analytical_eigenvalues, analytical_eigenvectors);
normalize_columns(analytical_eigenvectors);
sort_eigenpairs(analytical_eigenvalues, analytical_eigenvectors);
// Using self-made Jacobi solver
vec jacobi_eigenvalues;
mat jacobi_eigenvectors;
int iterations;
bool converged;
double eps = 1e-8;
int maxiter = 1e6;
jacobi_eigensolver(A, eps, jacobi_eigenvalues, jacobi_eigenvectors, maxiter, iterations, converged);
if (converged) {
cout << "Jacobi method converged in " << iterations << " iterations." << endl;
} else {
cout << "Jacobi method did not converge within the maximum number of iterations." << endl;
return 1;
}
normalize_columns(jacobi_eigenvectors); // Normalize the eigenvectors
sort_eigenpairs(jacobi_eigenvalues, jacobi_eigenvectors);
// Write results to file
ofstream outfile;
string filename = "eigen_data_" + to_string(n) + ".txt"; // e.g., eigen_data_10.txt
outfile.open(filename);
for(int i = 0; i < NUM_EIGENVECTORS; i++){
outfile << "analytical_solution_" << i+1 << "," << "jacobi_solution_" << i+1; // Column headers
if(i < NUM_EIGENVECTORS - 1){
outfile << ",";
} else {
outfile << endl;
}
}
write_boundary_conditions(outfile, NUM_EIGENVECTORS, DECIMAL_PRECISION); // Boundary condition at x=0
for(int j = 0; j < N; j++){
for(int i = 0; i < NUM_EIGENVECTORS; i++){
outfile << setprecision(DECIMAL_PRECISION) << scientific << analytical_eigenvectors(j,i) << "," << jacobi_eigenvectors(j,i);
if(i < NUM_EIGENVECTORS - 1){
outfile << ",";
} else {
outfile << endl;
}
}
}
write_boundary_conditions(outfile, NUM_EIGENVECTORS, DECIMAL_PRECISION); // Boundary condition at x=1
outfile.close();
return 0;
}
src/project2/generate_results.sh
Executable
+25
View File
@@ -0,0 +1,25 @@
#!/bin/bash
echo "Cleaning up old results..."
rm -f eigen_data_*.txt scaling_data*.txt
rm -f arma_eigenproblem_checker buckling_beam_solver off_diagonal_finder scaling_tester
rm -f ../../projects/project2/include/*.pdf
echo "Compiling the C++ code..."
g++ -larmadillo src/* arma_eigenproblem_checker.cpp -o arma_eigenproblem_checker
g++ -larmadillo src/* buckling_beam_solver.cpp -o buckling_beam_solver
g++ -larmadillo src/* off_diagonal_finder.cpp -o off_diagonal_finder
g++ -larmadillo src/* scaling_tester.cpp -o scaling_tester
echo "Compilation finished."
echo "Generating numerical results..."
./scaling_tester 250
./scaling_tester 250 dense
uv run python/scaling_plotter.py scaling_data.txt
uv run python/scaling_plotter.py scaling_data_dense.txt
./buckling_beam_solver 10
./buckling_beam_solver 100
uv run python/buckling_beam_plotter.py eigen_data_10.txt
uv run python/buckling_beam_plotter.py eigen_data_100.txt
echo "All done!"
src/project2/include/helpers.hpp
+18
View File
@@ -0,0 +1,18 @@
#include <armadillo>
#include <cmath>
arma::mat create_tridiagonal_matrix(int N, double h);
arma::mat create_dense_matrix(int N, double h);
void normalize_columns(arma::mat& A);
void analytical_solution(int N, double h, arma::vec& eigenvalues, arma::mat& eigenvectors);
double max_offdiag_symmetric(const arma::mat& A, int& k, int &l);
void jacobi_rotate(arma::mat& A, arma::mat& R, int k, int l);
void jacobi_eigensolver(const arma::mat& A, double eps, arma::vec& eigenvalues, arma::mat& eigenvectors, const int maxiter, int& iterations, bool& converged);
void sort_eigenpairs(arma::vec& eigenvalues, arma::mat& eigenvectors);
src/project2/off_diagonal_finder.cpp
+23
View File
@@ -0,0 +1,23 @@
#include "include/helpers.hpp"
#include <iostream>
#include <armadillo>
using namespace std;
using namespace arma;
int main(){
// Setup of the test matrix
mat A = {{1., 0., 0., 0.5},
{0., 1., -0.7, 0.},
{0., -0.7, 1., 0.},
{0.5, 0., 0., 1.}};
// Call the function
int k, l;
double max_val = max_offdiag_symmetric(A, k, l);
// Print the results
cout << "Max off-diagonal value: " << max_val << " at (" << k << ", " << l << ")" << endl;
return 0;
}
src/project2/python/buckling_beam_plotter.py
+38
View File
@@ -0,0 +1,38 @@
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import argparse
import os
parser = argparse.ArgumentParser(description='Plot buckling beam data.')
parser.add_argument('input_file', type=str, help='Path to the input CSV file')
args = parser.parse_args()
data = pd.read_csv(args.input_file)
x = np.linspace(0, 1, len(data))
fig, axs = plt.subplots(1, 2, figsize=(12, 5), sharex=True)
for col in data.columns:
c = f"C{int(col[-1])}"
ls = "-" if "jacobi" in col else "--"
label = col.replace("_", " ").title()
axs[0].plot(x, data[col], label=label, color=c, linestyle=ls)
if "analytical" in col:
jacobi_col = col.replace("analytical", "jacobi")
if jacobi_col in data.columns:
error = np.abs(data[col] - data[jacobi_col])
axs[1].plot(x, error, label=f'Error {label.split(" ")[-1]}', color=c, linestyle=':')
axs[1].set_xlabel('x')
axs[1].set_ylabel('Error')
axs[1].set_title('Buckling Beam Errors')
axs[1].legend()
axs[1].grid()
axs[0].set_xlabel('x')
axs[0].set_ylabel('Displacement')
axs[0].set_title('Buckling Beam Displacements')
axs[0].legend()
axs[0].grid()
dir = os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../projects/project2/include"))
filename = args.input_file.split("/")[-1].replace(".txt", "")
fig.savefig(f"{dir}/{filename}_plot.pdf")
src/project2/python/scaling_plotter.py
+37
View File
@@ -0,0 +1,37 @@
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import argparse
import os
parser = argparse.ArgumentParser(description='Plot scaling data.')
parser.add_argument('input_file', type=str, help='Path to the input CSV file')
parser.add_argument("--times", "-t", action="store_true", help="Plot times instead of iterations")
args = parser.parse_args()
data = pd.read_csv(args.input_file)
x = data["N"]
fig, ax = plt.subplots(figsize=(8, 5))
if args.times:
ax.plot(x, data["Jacobi_time"], label='Jacobi Time', color='C0', linestyle='-', marker='o')
ax.plot(x, data["Arma_time"], label='Armadillo Time', color='C1', linestyle='-', marker='o')
ax.set_ylabel('Time (s)')
ax.set_title('Scaling of Computation Time')
ax.set_yscale('log')
ax.legend()
ax.grid()
else:
ax.set_ylabel('Number of Iterations')
ax.set_title('Scaling of Jacobi Iterations')
ax.set_yscale('log')
ax.plot(x, data["Jacobi_iterations"], label='Jacobi Iterations', color='C0', linestyle='-', marker='o')
ax.grid()
ax.set_xlabel('Matrix Size N')
dir = os.path.abspath(os.path.join(os.path.dirname(__file__), "../../../projects/project2/include"))
filename = args.input_file.split("/")[-1].replace(".txt", "")
if args.times:
fig.savefig(f"{dir}/{filename}_times.pdf")
else:
fig.savefig(f"{dir}/{filename}_iterations.pdf")
src/project2/scaling_tester.cpp
+73
View File
@@ -0,0 +1,73 @@
#include <armadillo>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "include/helpers.hpp"
using namespace std;
using namespace arma;
int main(int argc, char** argv) {
int max_N = 200; // Maximum matrix size
if (argc > 1){ // Allow user to specify max matrix size during call
max_N = std::stoi(argv[1]);
}
string filename;
mat (*matrix_creator_func)(int, double);
if (argc > 2 && string(argv[2]) == "dense") { // Use dense matrix if specified
filename = "scaling_data_dense.txt";
matrix_creator_func = create_dense_matrix;
} else { // Default to tridiagonal matrix
filename = "scaling_data.txt";
matrix_creator_func = create_tridiagonal_matrix;
}
ofstream outfile;
outfile.open(filename);
outfile << "N,Arma_time,Analytical_time,Jacobi_time,Jacobi_iterations,Jacobi_converged\n";
for(int N = 10; N <= max_N; N += 10) { // Iterate over matrix sizes in steps of 10
double h = 1.0 / (N + 1);
mat A = matrix_creator_func(N, h);
// Armadillo solver
vec arma_eigenvalues;
mat arma_eigenvectors;
auto start_arma = chrono::high_resolution_clock::now();
eig_sym(arma_eigenvalues, arma_eigenvectors, A);
auto end_arma = chrono::high_resolution_clock::now();
chrono::duration<double> arma_duration = end_arma - start_arma;
// Analytical solver
vec analytical_eigenvalues;
mat analytical_eigenvectors;
auto start_analytical = chrono::high_resolution_clock::now();
analytical_solution(N, h, analytical_eigenvalues, analytical_eigenvectors);
auto end_analytical = chrono::high_resolution_clock::now();
chrono::duration<double> analytical_duration = end_analytical - start_analytical;
// Jacobi solver
vec jacobi_eigenvalues;
mat jacobi_eigenvectors;
int iterations;
bool converged;
double eps = 1e-8;
int maxiter = 1e6;
auto start_jacobi = chrono::high_resolution_clock::now();
jacobi_eigensolver(A, eps, jacobi_eigenvalues, jacobi_eigenvectors, maxiter, iterations, converged);
auto end_jacobi = chrono::high_resolution_clock::now();
chrono::duration<double> jacobi_duration = end_jacobi - start_jacobi;
// Write results to file
outfile << N << ","
<< arma_duration.count() << ","
<< analytical_duration.count() << ","
<< jacobi_duration.count() << ","
<< iterations << ","
<< (converged ? "true" : "false") << "\n";
cout << "Completed N = " << N << endl;
}
outfile.close();
return 0;
}
src/project2/src/helpers.cpp
+136
View File
@@ -0,0 +1,136 @@
#include "../include/helpers.hpp"
arma::mat create_tridiagonal_matrix(int N, double h) {
arma::mat A = arma::mat(N, N, arma::fill::zeros);
double diag_value = 2.0 / (h * h);
double off_diag_value = -1.0 / (h * h);
for (int i = 0; i < N; ++i) {
A(i, i) = diag_value;
if (i > 0) { // There is now lower diagnonal element in the first row
A(i, i - 1) = off_diag_value;
}
if (i < N - 1) { // There is no upper diagonal element in the last row
A(i, i + 1) = off_diag_value;
}
}
return A;
}
arma::mat create_dense_matrix(int N, double h) {
arma::mat A = arma::mat(N, N).randn();
A = arma::symmatu(A); // Make it symmetric
return A;
}
void normalize_columns(arma::mat& A) { // Useful for eigenvector normalization
A = arma::normalise(A, 2, 0); // 2-norm along columns (0)
// Flip sign to ensure first element is positive
for (size_t col = 0; col < A.n_cols; ++col) {
if (A(0, col) < 0) {
A.col(col) = -A.col(col);
}
}
}
void analytical_solution(int N, double h, arma::vec& eigenvalues, arma::mat& eigenvectors) {
eigenvalues = arma::vec(N);
eigenvectors = arma::mat(N, N);
double d = 2.0 / (h * h);
double a = -1.0 / (h * h);
for (int j = 0; j < N; ++j) {
eigenvalues(j) = d + 2 * a * std::cos((j + 1) * arma::datum::pi / (N + 1)); // Eigenvalues
for (int i = 0; i < N; ++i) {
eigenvectors(i, j) = std::sin((i + 1) * (j + 1) * arma::datum::pi / (N + 1)); // Eigenvectors
}
}
normalize_columns(eigenvectors); // Normalize the eigenvectors
}
double max_offdiag_symmetric(const arma::mat& A, int& k, int &l) {
int N = A.n_rows;
double max_val = 0.0;
for (int i = 0; i < N; ++i) {
for (int j = i + 1; j < N; ++j) { // Only check upper triangle
if (std::abs(A(i, j)) > max_val) {
max_val = std::abs(A(i, j));
k = i;
l = j;
}
}
}
return max_val;
}
void jacobi_rotate(arma::mat& A, arma::mat& R, int k, int l) {
if (A(k, l) != 0.0) {
double tau = (A(l, l) - A(k, k)) / (2.0 * A(k, l));
double t;
if (tau >= 0) {
t = 1.0 / (tau + std::sqrt(1.0 + tau * tau));
} else {
t = -1.0 / (-tau + std::sqrt(1.0 + tau * tau));
}
double c = 1.0 / std::sqrt(1 + t * t);
double s = t * c;
double temp;
// Update matrix A
double a_kk = A(k, k);
double a_ll = A(l, l);
A(k, k) = a_kk * c * c - 2.0 * A(k, l) * c * s + a_ll * s * s;
A(l, l) = a_ll * c * c + 2.0 * A(k, l) * c * s + a_kk * s * s;
A(k, l) = 0.0; // Hard-coding to zero
A(l, k) = 0.0; // Hard-coding to zero
for (int i = 0; i < A.n_rows; ++i) {
if (i != k && i != l) {
temp = A(i, k);
A(i, k) = temp * c - A(i, l) * s;
A(k, i) = A(i, k); // Since A is symmetric
A(i, l) = temp * s + A(i, l) * c;
A(l, i) = A(i, l); // Since A is symmetric
}
// Update the eigenvector matrix R
temp = R(i, k);
R(i, k) = temp * c - R(i, l) * s;
R(i, l) = temp * s + R(i, l) * c;
}
}
}
void jacobi_eigensolver(const arma::mat& A, double eps, arma::vec& eigenvalues, arma::mat& eigenvectors,
const int maxiter, int& iterations, bool& converged) {
int N = A.n_rows;
arma::mat A_copy = A; // Make a copy to avoid modifying the original matrix
eigenvectors = arma::eye<arma::mat>(N, N); // Initialize eigenvector matrix as identity
iterations = 0;
converged = false;
int k, l;
double max_offdiag;
while (iterations < maxiter) {
max_offdiag = max_offdiag_symmetric(A_copy, k, l);
if (max_offdiag < eps) {
converged = true;
break;
}
jacobi_rotate(A_copy, eigenvectors, k, l);
iterations++;
}
eigenvalues = A_copy.diag(); // Extract eigenvalues from the diagonal of A_copy
}
void sort_eigenpairs(arma::vec& eigenvalues, arma::mat& eigenvectors) {
arma::uvec indices = arma::sort_index(eigenvalues, "ascend");
eigenvalues = eigenvalues(indices);
eigenvectors = eigenvectors.cols(indices);
}