-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathutils.h
111 lines (86 loc) · 2.46 KB
/
utils.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
#pragma once
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
void afficher(vector <double>&v) {
cout << "debut de l'affichage" << endl;
for (int i = 0; i < v.size(); i = i + 1) {
cout << v[i] << " ";
};
cout << endl;
}
double fRand(double fMin, double fMax)
{
double f = (double)rand() / RAND_MAX;
return fMin + f * (fMax - fMin);
}
void print(vector<int> const& v)
{
for (auto& i : v)
std::cout << i << ' ';
std::cout << '\n';
}
void firststepmodif(vector <double> &a, vector <double>& b, vector <double>& c, vector <double>& f ) {// premiere etape de la modification de la matrix
vector<double> temp;
temp.resize(4);
for (int j = 1; j < a.size(); j++) {
temp.clear();
temp.resize(4);
double zeroavoid;
if (b[j - 1] == 0) {//pour eviter de diviser pa zero on divisera par 1 en considerant qu'il s'agit d'une operation matricielle
//et que cette division ne se ferai pas si le coeficient bJ est deja nul
zeroavoid = 1;
}
else {
zeroavoid = b[j - 1];
}
temp[0] = a[j - 1] * a[j]/ zeroavoid;
temp[1] = b[j] - c[j - 1] * a[j] / zeroavoid;
temp[2] = c[j];
temp[3] = f[j] - f[j - 1] * a[j]/ zeroavoid;
a[j] = temp[0];
b[j] = temp[1];
c[j] = temp[2];
f[j] = temp[3];
};
}
void secondstepmodif(vector <double>& a, vector <double>& b, vector <double>& c, vector <double>& f) { // deuxieme etape de la modification de la matrix
vector<int> temp;
temp.resize(4);
for (int j = a.size()-3; j >=0; j--) {
temp.clear();
temp.resize(4);
double zeroavoid;
if (b[j +1] == 0) {//pour eviter de diviser pa zero on divisera par 1 en considerant qu'il s'agit d'une operation matricielle
//et que cette division ne se ferai pas si le coeficient bJ est deja nul
zeroavoid = 1;
}
else {
zeroavoid = b[j + 1];
}
temp[0] = a[j] -a[j+1] * c[j] / zeroavoid;
temp[1] = b[j];
temp[2] = c[j]*c[j+1]/ zeroavoid;
temp[3] = f[j] - f[j + 1] * c[j] / zeroavoid;
a[j] = temp[0];
b[j] = temp[1];
c[j] = temp[2];
f[j] = temp[3];
}
}
vector<double>tridiagonal_solver(vector<double>& a, vector<double>& b, vector<double>& c, vector<double>& f) {
int n = f.size();
vector<double> x(n);
for (int i = 1; i < n; i++) {
double m = a[i - 1] / b[i - 1];
b[i] -= m * c[i - 1];
f[i] -= m * f[i - 1];
}
// solve for last x value
x[n - 1] = f[n - 1] / b[n - 1];
// solve for remaining x values by back substitution
for (int i = n - 2; i >= 0; i--)
x[i] = (f[i] - c[i] * x[i + 1]) / b[i];
return x;
}