Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <vector>
- #include <math.h>
- #include <complex>
- #include <stdlib.h>
- #include <sstream>
- #include <iomanip>
- using namespace std;
- class Vector;
- class Matrix {
- protected:
- vector<vector<double>> data;
- int n;
- int m;
- public:
- Matrix() {
- n = 0;
- m = 0;
- }
- Matrix (int x, int y) {
- n = x;
- m = y;
- data.assign(n, vector<double>(m, 0));
- }
- Matrix& operator+ (const Matrix& matrix) {
- if (matrix.m != m || matrix.n != n) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] = data[i][j] + matrix.data[i][j];
- }
- }
- return *ans;
- }
- Matrix& operator- (const Matrix& matrix) {
- if (matrix.m != m || matrix.n != n) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] = data[i][j] - matrix.data[i][j];
- }
- }
- return *ans;
- }
- Matrix& operator* (const Matrix& matrix) {
- if (matrix.n != m) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, matrix.m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- for (int k = 0; k < m; k++) {
- ans->data[i][j] += data[i][k]*matrix.data[k][j];
- }
- }
- }
- return *ans;
- }
- Matrix& operator* (double a) {
- Matrix *ans = new Matrix(n, m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] += data[i][j]*a;
- }
- }
- return *ans;
- }
- Matrix& transpose() {
- Matrix *ans = new Matrix(m, n);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[j][i] = data[i][j];
- }
- }
- return *ans;
- }
- Matrix& operator| (const Matrix& matrix) {
- if (matrix.n != n) {
- throw runtime_error("Error: the dimensional problem occurred");
- }
- Matrix *ans = new Matrix(n, matrix.m + m);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- ans->data[i][j] = data[i][j];
- }
- for (int j = m; j < matrix.m + m; j++) {
- ans->data[i][j] = matrix.data[i][j-m];
- }
- }
- return *ans;
- }
- vector<double>& operator[] (int x) {
- if (x >= n) {
- throw runtime_error("Error: index out of range");
- }
- return data[x];
- }
- friend istream& operator>> (istream &is, Matrix &matrix);
- friend ostream& operator<< (ostream &os, Matrix &matrix);
- double get(int i, int j) {
- if (i >= n || j >= m) {
- throw out_of_range("Out of bound for get");
- }
- return data[i][j];
- }
- bool is_square() {
- return n == m;
- }
- int get_n() {
- return n;
- }
- int get_m() {
- return m;
- }
- void normalize_diagonally() {
- for (int i = 0; i < n; i++) {
- if (data[i][i] != 0) {
- for (int j = i + 1; j < m; j++) {
- data[i][j] /= data[i][i];
- }
- data[i][i] = 1;
- }
- }
- }
- void separate() {
- // for (int i = 0; i < n; i++) {
- // for (int j = 0; j < m - 1; j++) {
- // // if (abs(data[i][j]) < 0.005) {
- // // // // cout << fixed << abs(data[i][j]) << ' ';
- // // } else
- // // // // cout << fixed << data[i][j] << ' ';
- // }
- // // cout << '\n';
- // }
- // for (int j = 0; j < n; j++) {
- // if (abs(data[j][m - 1]) < 0.005) {
- // // cout << fixed << abs(data[j][m - 1]) << ' ';
- // } else
- // // cout << fixed << data[j][m - 1] << ' ';
- // }
- // // cout << '\n';
- }
- void last() {
- // for (int j = 0; j < n; j++) {
- // if (abs(data[j][m - 1]) < 0.005) {
- // // cout << fixed << abs(data[j][m - 1]) << ' ';
- // } else
- // // cout << fixed << data[j][m - 1] << ' ';
- // }
- // // cout << '\n';
- }
- Matrix& inverse();
- operator Vector();
- };
- class Vector : public Matrix {
- public:
- double normalize() {
- double l = 0;
- for (int i = 0; i < n; i++) {
- l += data[i][0] * data[i][0];
- }
- return l;
- }
- Vector(vector<double> input) : Matrix(input.size(), 1) {
- for (int i = 0 ; i < input.size(); i++) {
- data[i][0] = input[i];
- }
- }
- Vector() {
- n = 0;
- m = 0;
- }
- double operator[] (int i) {
- if (i >= n) {
- throw runtime_error("Error: Out of range");
- }
- return data[i][0];
- }
- friend istream& operator>> (istream &is, Vector& vector);
- };
- Matrix::operator Vector() {
- if (m != 1) {
- throw runtime_error("Error in cast from Matrix to Vector");
- }
- vector<double> a;
- for (auto el : data) {
- a.push_back(el[0]);
- }
- return *(new Vector(a));
- }
- class LeastSquareAproxMatrix : public Matrix {
- public:
- LeastSquareAproxMatrix(int m, vector<double> input) : Matrix(input.size(), m + 1) {
- for (int i = 0; i < input.size(); i++) {
- double pow = 1;
- for (int j = 0; j <= m; j++) {
- data[i][j] = pow;
- pow *= input[i];
- }
- }
- }
- };
- class SquareMatrix : public Matrix {
- public:
- SquareMatrix() {
- n = 0;
- m = 0;
- }
- SquareMatrix* upper_triangular() {
- SquareMatrix* ans = new SquareMatrix(n);
- for (int i = 0; i < n; i++) {
- for (int j = i + 1; j < n; j++) {
- (*ans)[i][j] = data[i][j];
- }
- }
- return ans;
- }
- SquareMatrix* lower_triangular() {
- SquareMatrix* ans = new SquareMatrix(n);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < i; j++) {
- (*ans)[i][j] = data[i][j];
- }
- }
- return ans;
- }
- SquareMatrix& build_apha() {
- SquareMatrix* ans = new SquareMatrix(n);
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- if (i == j) {
- (*ans)[i][j] = 0;
- } else {
- (*ans)[i][j] = -data[i][j]/data[i][i];
- }
- }
- }
- return *ans;
- }
- Vector& build_beta(Vector& b) {
- vector<double> a;
- for (int i = 0; i < n; i++) {
- a.push_back(b[i] / data[i][i]);
- }
- Vector* ans = new Vector(a);
- return *ans;
- }
- bool is_diagonal_dominant() {
- for (int i = 0; i < n; i++) {
- double sum = 0;
- for (int j = 0; j < m; j++) {
- if (i != j) {
- sum += data[i][j];
- }
- }
- if (data[i][i] < sum) {
- return 0;
- }
- }
- return 1;
- }
- SquareMatrix(int x) : Matrix(x, x) {};
- friend istream& operator>> (istream &is, SquareMatrix &matrix);
- void modify_for_l() {
- for (int i = 0; i < n; i++) {
- for (int j = 0; j <= i; j++) {
- data[i][j] = 1;
- }
- }
- }
- };
- class IdentityMatrix : public SquareMatrix {
- public:
- IdentityMatrix() {
- n = 0;
- m = 0;
- }
- IdentityMatrix(int x) : SquareMatrix(x) {
- for (int i = 0; i < x; i++) {
- data[i][i] = 1;
- }
- }
- };
- class EliminationMatrix : public IdentityMatrix {
- public:
- EliminationMatrix(Matrix& matrix, int j, int i) : IdentityMatrix(matrix.get_n()) {
- data[j - 1][i - 1] = -matrix.get(j - 1, i - 1) / matrix.get(i - 1, i - 1);
- }
- void negative() {
- for(int i = 0; i < n; i++) {
- for (int j = 0; j < m; j++) {
- data[i][j] = -data[i][j];
- }
- }
- }
- };
- class PermutationMatrix : public IdentityMatrix {
- public:
- PermutationMatrix(int x, int i, int j) : IdentityMatrix(x) {
- data[j - 1][j - 1] = 0;
- data[i - 1][i - 1] = 0;
- data[i - 1][j - 1] = 1;
- data[j - 1][i - 1] = 1;
- }
- };
- istream& operator>> (istream &is, SquareMatrix& matrix) {
- is >> matrix.n;
- matrix.m = matrix.n;
- matrix.data.assign(matrix.n, vector<double>(matrix.m));
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- is >> matrix.data[i][j];
- }
- }
- return is;
- }
- istream& operator>> (istream &is, Vector& matrix) {
- is >> matrix.n;
- matrix.m = 1;
- matrix.data.assign(matrix.n, vector<double>(matrix.m));
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- is >> matrix.data[i][j];
- }
- }
- return is;
- }
- istream& operator>> (istream &is, Matrix& matrix) {
- is >> matrix.n >> matrix.m;
- matrix.data.assign(matrix.n, vector<double>(matrix.m));
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- is >> matrix.data[i][j];
- }
- }
- return is;
- }
- ostream& operator<< (ostream &os, Matrix& matrix) {
- for (int i = 0; i < matrix.n; i++) {
- for (int j = 0; j < matrix.m; j++) {
- if (abs(matrix.data[i][j]) < 0.005) {
- os << fixed << abs(matrix.data[i][j]) << ' ';
- } else
- os << fixed << matrix.data[i][j] << ' ';
- }
- os << '\n';
- }
- return os;
- }
- bool gaus_elimination(Matrix& a, int &step) {
- bool sign = 1;
- for (int i = 0; i < a.get_n(); i++) {
- double m = 0;
- int m_i = 0;
- for (int j = i; j < a.get_n(); j++) {
- if (abs(a.get(j, i)) > abs(m)) {
- m = a.get(j, i);
- m_i = j;
- }
- }
- if (m_i != i) {
- // // cout << fixed << "step #" << step << ": permutation" << '\n';
- PermutationMatrix p(a.get_n(), m_i + 1, i + 1);
- a = p * a;
- // a.separate();
- sign = !sign;
- step++;
- }
- for (int j = i + 1; j < a.get_n(); j++) {
- if (a.get(j, i) != 0) {
- // // cout << fixed << "step #" << step << ": elimination" << '\n';
- step++;
- EliminationMatrix e(a, j + 1, i + 1);
- a = e * a;
- // a.separate();
- }
- }
- }
- return sign;
- }
- void gaus_back_elimination(Matrix& a, int &step) {
- for (int i = a.get_n() - 1; i >= 0; i--) {
- if (a.get(i, i) != 0) {
- for (int j = i - 1 ; j >= 0; j--) {
- if (a.get(j, i) != 0) {
- // // cout << fixed << "step #" << step << ": elimination" << '\n';
- step++;
- EliminationMatrix e(a, j + 1, i + 1);
- a = e * a;
- a.separate();
- }
- }
- }
- }
- }
- Matrix& Matrix::inverse() {
- Matrix* ans = new Matrix(n, m);
- IdentityMatrix I(n);
- Matrix aug = *this | I;
- int step = 1;
- gaus_elimination(aug, step);
- gaus_back_elimination(aug, step);
- aug.normalize_diagonally();
- for (int j = 0; j < aug.get_n(); j++) {
- for (int i = aug.get_n(); i < aug.get_m(); i++) {
- (*ans)[j][i - aug.get_n()] = aug.get(j, i);
- }
- }
- return *ans;
- }
- class ColumnVector : public Matrix {
- public:
- ColumnVector() {
- n = 0;
- m = 0;
- }
- ColumnVector(int x) {
- n = x;
- m = 1;
- data.resize(n, vector<double>(1));
- }
- friend istream& operator>> (istream &is, ColumnVector &matrix);
- friend ostream& operator<< (ostream &os, ColumnVector &matrix);
- };
- istream& operator>> (istream &is, ColumnVector& matrix) {
- is >> matrix.n;
- matrix.m = 1;
- matrix.data.resize(matrix.n,vector<double>(1));
- for (int i = 0; i < matrix.n; i++) {
- is >> matrix.data[i][0];
- }
- return is;
- }
- ostream& operator<< (ostream &os, ColumnVector& matrix) {
- for (int i = 0; i < matrix.n; i++) {
- os << matrix.data[i][0] << ' ';
- }
- os << '\n';
- return os;
- }
- bool is_singular( Matrix a) {
- for (int i = 0; i < a.get_n(); i++) {
- double m = 0;
- int m_i = 0;
- for (int j = i; j < a.get_n(); j++) {
- if (abs(a.get(j, i)) > abs(m)) {
- m = a.get(j, i);
- m_i = j;
- }
- }
- if (m_i != i) {
- PermutationMatrix p(a.get_n(), m_i + 1, i + 1);
- a = p * a;
- }
- for (int j = i + 1; j < a.get_n(); j++) {
- if (a.get(j, i) != 0) {
- EliminationMatrix e(a, j + 1, i + 1);
- a = e * a;
- }
- }
- }
- for (int i = 0; i < a.get_n(); i++) {
- if (a.get(i, i) == 0) {
- return 0;
- }
- }
- return 1;
- }
- double calc_epsilon(Vector& a, Vector& b) {
- if (a.get_n() != b.get_n()) {
- throw runtime_error("cannot calculate epsilon");
- }
- double ans = 0;
- for (int i = 0; i < a.get_n(); i++) {
- ans += (a[i] - b[i]) * (a[i] - b[i]);
- }
- return sqrt(ans);
- }
- void jacobi_method(SquareMatrix& alpha, Vector& beta, Vector& x, double epsilon, int step) {
- cout << "x(" << step << "):\n";
- Vector new_x = alpha*x + beta;
- cout << new_x;
- cout << "e:\n" << calc_epsilon(new_x, x) << '\n';
- if (calc_epsilon(new_x, x) < epsilon) {
- x = new_x;
- return;
- }
- jacobi_method(alpha, beta, new_x, epsilon, step + 1);
- x = new_x;
- }
- double min(vector<double>& a) {
- double m = 1000000000;
- for (double el : a) {
- m = min(el, m);
- }
- return m;
- }
- double max(vector<double>& a) {
- double m = -100000000;
- for (double el : a) {
- m = max(el, m);
- }
- return m;
- }
- void seidel_method(Matrix& B, Matrix& C, Vector& beta, Vector& x, double epsilon, int step) {
- cout << "x(" << step << "):\n";
- IdentityMatrix I(B.get_n());
- Vector new_x = (I - B).inverse() * C * x + (I - B).inverse() * beta;
- cout << new_x;
- cout << "e:\n" << calc_epsilon(new_x, x) << '\n';
- if (calc_epsilon(new_x, x) < epsilon) {
- x = new_x;
- return;
- }
- seidel_method(B, C, beta, new_x, epsilon, step + 1);
- x = new_x;
- }
- FILE* pipe;
- int main() {
- pipe = popen("/usr/bin/gnuplot -persist", "w");
- cout.precision(2);
- double v, k;
- cin >> v;
- cin >> k;
- double a1, b1, a2, b2;
- cin >> a1 >> b1 >> a2 >> b2;
- double T;
- cin >> T;
- double n;
- cin >> n;
- SquareMatrix A(2);
- A[0][1] = -a2*b1/b2;
- A[1][0] = a1*b2/b1;
- vector<double> t;
- vector<double> ks;
- vector<double> vs;
- Vector X_0({v, k});
- // real solution
- // for (double i = 0; i <= T; i += T / n) {
- // t.push_back(i);
- // complex<double> l1 = sqrt(abs(A[0][1] * A[1][0]))*1i;
- // complex<double> l2 = -sqrt(abs(A[0][1] * A[1][0]))*1i;
- // complex<double> y = (pow(M_E, l2*i) - pow(M_E, l1*i)) / (l2 - l1);
- // complex<double> x = pow(M_E, l1*i)-l1*y;
- // IdentityMatrix I(2);
- // Matrix V = I*x.real() + A*y.real();
- // Matrix Ans = V * X_0;
- // vs.push_back(Ans[0][0]);
- // ks.push_back(Ans[1][0]);
- // }
- // just formula
- for (double i = 0; i <= T; i += T / n) {
- vs.push_back(v*cos(sqrt(a1 * a2) * i) - k*(sqrt(a2)*b1/b2*sqrt(a1))*sin(sqrt(a1*a2)*i));
- ks.push_back(v*(sqrt(a1)*b2/b1*sqrt(a2))*sin(sqrt(a1*a2)*i)+k*cos(sqrt(a1*a2)*i));
- t.push_back(i);
- }
- cout << "t:\n";
- for (double el : t) {
- cout << fixed << el << ' ';
- }
- fputs(("plot ["+ to_string(0 - 1) + ":" + to_string(int (T + 1)) + "]["+ to_string(int(min(min(ks), min(vs)) - 1)) + ":" + to_string(int(max(max(ks), max(vs))+ 1)) + "] '-' using 1:2 title 'v(t)' with lines , '-' title 'k(t)' with lines").c_str(), pipe);
- fprintf(pipe, "\n");
- cout << "\nv:\n";
- for (int i = 0; i < vs.size(); i++) {
- cout << fixed << vs[i] << ' ';
- fprintf(pipe, "%f\t%f\n", t[i], vs[i]);
- }
- fprintf(pipe, "%s\n", "e");
- cout << "\nk:\n";
- for (int i = 0; i < vs.size(); i++) {
- cout << fixed << ks[i] << ' ';
- fprintf(pipe, "%f\t%f\n", t[i], ks[i]);
- }
- fprintf(pipe, "%s\n", "e");
- fputs(("plot ["+ to_string(int(min(ks) - 1)) + ":" + to_string(int(max(ks) + 1)) + "]["+ to_string(int(min(vs) - 1)) + ":" + to_string(int(max(vs) + 1)) + "] '-' using 1:2 title 'v(k)' with lines").c_str(), pipe);
- fprintf(pipe, "\n");
- for (int i = 0; i < vs.size(); i++) {
- cout << fixed << ks[i] << ' ';
- fprintf(pipe, "%f\t%f\n", ks[i], vs[i]);
- }
- fprintf(pipe, "%s\n", "e");
- pclose(pipe);
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement