Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <cmath>
- using namespace std;
- double matrix[10010][10010];
- void buildMatrixForFirst(int n, int m, double hY, double hX)
- {
- for(int i = 0; i<n;++i)
- for(int j=0;j<m;++j)
- {
- matrix[i * m + j][i * m + j] = 2 / hY / hY;
- if(i > 0)
- {
- matrix[i * m + j][(i - 1) * m + j] = -1 / hX / hX - 1 / hY / hY ;
- }
- if(i < n - 1)
- {
- matrix[i * m + j][(i + 1) * m + j] = -1 / hX / hX - 1 / hY / hY ;
- }
- if(j > 0 && i > 0)
- {
- matrix[i * m + j][(i - 1) * m + (j - 1)] = -1 / hX / hX / 2;
- }
- if(j > 0 && i < n - 1)
- {
- matrix[i * m + j][(i + 1) * m + (j - 1)] = -1 / hX / hX / 2;
- }
- if(j < m - 1 && i > 0)
- {
- matrix[i * m + j][(i - 1) * m + (j + 1)] = -1 / hX / hX / 2;
- }
- if(j < m - 1 && i < n - 1)
- {
- matrix[i * m + j][(i + 1) * m + (j + 1)] = -1 / hX / hX / 2;
- }
- matrix[i * m + j][n * m] = 1;
- }
- }
- void buildMatrixForSecond(int n, int m, double hY, double hX)
- {
- for(int i = 0; i<n;++i)
- for(int j=0;j<m;++j)
- {
- matrix[i * m + j][i * m + j] = 2 / hY / hY;
- if(i > 0)
- {
- matrix[i * m + j][(i - 1) * m + j] = -1 / hX / hX - 1 / hY / hY ;
- }
- if(i < n - 1)
- {
- matrix[i * m + j][(i + 1) * m + j] = -1 / hX / hX - 1 / hY / hY ;
- }
- if(j > 0 && i > 0)
- {
- matrix[i * m + j][(i - 1) * m + (j - 1)] = -1 / hX / hX;
- }
- if(j < m - 1 && i > 0)
- {
- matrix[i * m + j][(i - 1) * m + (j + 1)] = -1 / hX / hX;
- }
- matrix[i * m + j][n * m] = 1;
- }
- }
- int now[10000];
- bool solve(int n)
- {
- for(int i = 0; i<n;++i)
- {
- if(abs(matrix[i][i]) < 1e-9)
- {
- //cout << "LOOOOOOL";
- bool find = 0;
- for(int j = i + 1; j < n; ++j)
- {
- if(abs(matrix[j][i]) > 1e-9)
- {
- find = 1;
- for(int u = 0; u < n + 1; ++u)
- swap(matrix[i][u], matrix[j][u]);
- swap(now[i], now[j]);
- break;
- }
- }
- if(!find)
- return 0;
- }
- double div = matrix[i][i];
- for(int j=0;j<n + 1;++j)
- matrix[i][j] /= div;
- for(int j=0;j<n;++j)
- {
- if(j != i)
- {
- double mnoj = matrix[j][i];
- for(int u = 0; u < n + 1; ++u)
- {
- matrix[j][u] -= matrix[i][u] * mnoj;
- }
- }
- }
- }
- return 1;
- }
- double u[10010][10010];
- double newU[10010][10010];
- double nextTurn(int a, int b, double h)
- {
- double eps = 0;
- for (int i = 2; i < a; ++i)
- for (int j = 2; j < b; ++j)
- {
- newU[i][j] = (u[i-1][j] + u[i][j - 1] + u[i+1][j] + u[i][j+1] + h*h) / 4;
- eps = max(eps, abs(newU[i][j] - u[i][j]));
- }
- swap(u, newU);
- return eps;
- }
- int main()
- {
- double a, b, hX, hY, eps;
- cin>>a>>b;
- cin >> hX >> hY;
- eps = 1e-9;
- cout.precision(16);
- int n = a/hX, m = b/hY;
- buildMatrixForFirst(n - 2, m - 2, hX, hY);
- //buildMatrixForSecond(n - 2, m - 2, hX, hY);
- for (int i = 0; i < 10000; ++i) {
- now[i] = i;
- }
- if(solve((n - 2) * (m - 2)))
- cout<<"KEK";
- else
- cout<<"KEEEEEEEEK";
- for(int i = 0; i < n - 2; ++i)
- {
- for(int j = 0; j < m - 2; ++j)
- cout<<matrix[now[i * (m - 2) + j]][(n - 2) * (m - 2)]<<" ";
- cout<<endl;
- }
- cout<<endl;
- //while(nextTurn(n, m, h) > eps);
- //for(int i = 2; i < n; ++i)
- //{
- // for(int j = 2; j < m; ++j)
- // cout<<u[i][j]<<" ";
- // cout<<endl;
- //}
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement