sim_ressorts/labo_physique/Solvers.h

148 lines
3.7 KiB
C
Raw Normal View History

2024-02-27 13:20:47 -05:00
#pragma once
/**
* @file Solvers.hpp
*
* @brief Implémentation de plusieurs algorihtmes de solveurs pour un système
* d'équations linéaires
*
* Nom:
* Code permanent :
* Email :
*
*/
#include "Math3D.h"
2024-03-09 16:19:14 -05:00
namespace gti320 {
2024-02-27 13:20:47 -05:00
// Identification des solveurs
2024-03-09 16:19:14 -05:00
enum eSolverType {
kNone, kGaussSeidel, kColorGaussSeidel, kCholesky
};
2024-02-27 13:20:47 -05:00
// Paramètres de convergences pour les algorithmes itératifs
static const float eps = 1e-4f;
static const float tau = 1e-5f;
/**
* Résout Ax = b avec la méthode Gauss-Seidel
*/
2024-03-09 16:19:14 -05:00
static void gaussSeidel(const Matrix<float, Dynamic, Dynamic> &A,
const Vector<float, Dynamic> &b,
Vector<float, Dynamic> &x, int k_max) {
2024-02-27 13:20:47 -05:00
// TODO
//
// Implémenter la méthode de Gauss-Seidel
2024-03-09 16:19:14 -05:00
int n = b.size();
x.resize(n);
x.setZero();
bool converged = false;
int k = 0;
do {
Vector<float, Dynamic> nx = x;
for (int i = 0; i < n; i++) {
nx(i) = b(i);
for (int j = 0; j < i; j++) {
nx(i) = nx(i) - A(i, j) * nx(j);
}
for (int j = i + 1; j < n; j++) {
nx(i) = nx(i) - A(i, j) * x(j);
}
nx(i) = nx(i) / A(i, i);
}
k++;
Vector<float, Dynamic> r = A * x - b;
converged = k >= k_max ||
(nx - x).norm() / nx.norm() < tau ||
r.norm() / b.norm() < eps;
x = nx;
} while (!converged);
2024-02-27 13:20:47 -05:00
}
/**
* Résout Ax = b avec la méthode Gauss-Seidel (coloration de graphe)
*/
2024-03-09 16:19:14 -05:00
static void gaussSeidelColor(const Matrix<float, Dynamic, Dynamic> &A, const Vector<float, Dynamic> &b,
Vector<float, Dynamic> &x, const Partitions &P, const int maxIter) {
2024-02-27 13:20:47 -05:00
// TODO
//
// Implémenter la méthode de Gauss-Seidel avec coloration de graphe.
// Les partitions avec l'index de chaque particule sont stockées dans la table des tables, P.
}
/**
* Résout Ax = b avec la méthode de Cholesky
*/
2024-03-09 16:19:14 -05:00
static void cholesky(const Matrix<float, Dynamic, Dynamic> &A,
const Vector<float, Dynamic> &b,
Vector<float, Dynamic> &x) {
int n = A.rows();
x.resize(n);
x.setZero();
2024-02-27 13:20:47 -05:00
// TODO
//
// Calculer la matrice L de la factorisation de Cholesky
2024-03-09 16:19:14 -05:00
auto L = Matrix<float, Dynamic, Dynamic>(n, n);
for (int j = 0; j < n; j++) {
for (int i = j; i < n; i++) {
float s = 0;
2024-02-27 13:20:47 -05:00
2024-03-09 16:19:14 -05:00
for (int k = 0; k < j; k++) {
s += L(i, k) * L(j, k);
}
2024-02-27 13:20:47 -05:00
2024-03-09 16:19:14 -05:00
if (i == j) {
L(i, i) = std::sqrt(A(i, i) - s);
} else {
L(i, j) = (A(i, j) - s) / L(j, j);
}
}
}
2024-02-27 13:20:47 -05:00
// TODO
//
// Résoudre Ly = b
2024-03-09 16:19:14 -05:00
auto y = Vector<float, Dynamic>(n);
for (int i = 0; i < n; i++) {
y(i) = b(i);
2024-02-27 13:20:47 -05:00
2024-03-09 16:19:14 -05:00
for (int j = 0; j < i; j++) {
y(i) -= L(i, j) * y(j);
}
2024-02-27 13:20:47 -05:00
2024-03-09 16:19:14 -05:00
y(i) /= L(i, i);
}
2024-02-27 13:20:47 -05:00
// TODO
//
// Résoudre L^t x = y
//
2024-03-09 16:19:14 -05:00
// Remarque : ne pas calculer la transposée de L, c'est inutilement
2024-02-27 13:20:47 -05:00
// coûteux.
2024-03-09 16:19:14 -05:00
for (int i = n - 1; i >= 0; i--) {
x(i) = y(i);
2024-02-27 13:20:47 -05:00
2024-03-09 16:19:14 -05:00
for (int j = i + 1; j < n; j++) {
x(i) -= L(j, i) * x(j);
}
2024-02-27 13:20:47 -05:00
2024-03-09 16:19:14 -05:00
x(i) /= L(i, i);
}
2024-02-27 13:20:47 -05:00
}
}