#include using namespace std; #include #include #include #include #include #include #define PI 3.14159265358979323846 using namespace std; void fill_matrix(float **matrix, int n, float t, float h, float *vector_u, float *vector_v, float *vector_b) { for (int i = 0; i < n; i ++) { for (int j = 0; j < n; j ++) { matrix[i][j] = 0; } } for (int i = 0; i < n; i ++) { //столбец b if (i < n/2) { vector_b[i] = vector_v[i]/t; } else { vector_b[i] = vector_u[i-(n/2)+1]/t; } for (int j = 0; j < n; j ++) { //диагональ if (i == j) { matrix[i][j] = 1/t; } if (i > 0 && i < (n/2)-1) { matrix[i][i-1] = -vector_u[i]/h; matrix[i][i+1] = vector_u[i]/h; } if (i == (n/2)-1) { matrix[i][i-1] = -vector_u[i]/h; } if (i == n/2) { matrix[i][i+1] = vector_v[i-(n/2)+1]/h; } if (i > (n/2) && i < n-1) { matrix[i][i-1] = -vector_v[i-(n/2)+1]/h; matrix[i][i+1] = vector_v[i-(n/2)+1]/h; } //поддиагональная диагональ :D if (i >= n/2 && i < n-2) { matrix[i][i-(n/2)] = -(vector_u[i-(n/2)+1]+vector_v[i-(n/2)+1])/h; matrix[i][i-(n/2)+2] = (vector_u[i-(n/2)+1]+vector_v[i-(n/2)+1])/h; } if (i == n-2) { matrix[i][i-(n/2)] = -(vector_u[i-(n/2)+1]+vector_v[i-(n/2)+1])/h; } if (i == n-1) { matrix[i][i-(n/2)-1] = vector_u[n/2]/h; matrix[i][i-(n/2)] = -4*vector_u[n/2]/h; } //наддиагональная диагональ :D if (i == 0) { matrix[i][i+(n/2)] = 4*vector_v[i]/h; matrix[i][i+(n/2)+1] = -vector_v[i]/h; } if (i == 1) { matrix[i][i+(n/2)] = (vector_u[i]+vector_v[i])/h; } if (i > 1 && i < (n/2)) { matrix[i][i+(n/2)-2] = -(vector_u[i]+vector_v[i])/h; matrix[i][i+(n/2)] = (vector_u[i]+vector_v[i])/h; } } } return; } std::ofstream ofs ("output.txt", std::ofstream::out); void print_matrix(float **matrix, int n) { for (int i = 0; i < n; i ++) { for (int j = 0; j < n; j ++) { ofs << matrix[i][j] << setw(13) << ' ' ; } ofs << endl; } ofs << endl; return; } void print_vector(float *p_vector, int n) { for (int i = 0; i < n; i ++) { ofs << p_vector[i] << ' ' ; } ofs << endl; } float** matrix_return(float cos, float sin, int n, int c, int u) { int i, j; float** mat = (float**)malloc(n * sizeof(float*)); for (i = 0; i < n; i ++) { mat[i] = (float*)malloc(n * sizeof(float)); } for (i = 0; i < n; i ++) { for (j = 0; j < n; j ++) { mat[i][j] = 0; } mat[i][i] = 1; } mat[c][c] = cos; mat[c][u+1] = -sin; mat[u+1][c] = sin; mat[u+1][u+1] = cos; return (mat); } float** multiplication(float **mat1, float **mat2, int n) { int i, j, k; float** mat3 = (float**)malloc(n * sizeof(float*)); for (i = 0; i < n; i ++) { mat3[i] = (float*)malloc(n * sizeof(float)); } for (i = 0; i < n; i ++) { for (j = 0; j < n; j ++) { mat3[i][j] = 0; } } for (i = 0; i < n; i ++) { for (j = 0; j < n; j ++) { for (k = 0; k < n; k ++) { mat3[i][j] += mat1[i][k] * mat2[k][j]; } } } return (mat3); } float* multiplication1(float **mat1, float *m_vector, int n) { int i, j, k; float* mat3 = (float*)malloc(n*sizeof(float)); for (i = 0; i < n; i ++) { mat3[i] = 0; } for (i = 0; i < n; i ++) { for (k = 0; k < n; k ++) { mat3[i] += mat1[i][k] * m_vector[k]; } } return(mat3); } float** fill_matrix_for_reverse(float **buf, float **single, float **res, int n, int c) { int i, j, k; float t, sum; float** buf1 = (float**)malloc(n * sizeof(float*)); for (i = 0; i < n; i ++) { buf1[i] = (float*)malloc(n * sizeof(float)); } for (i = 0; i < n; i ++) { for (j = 0; j < n; j++) { buf1[i][j] = buf[i][j]; } } for (j = n - 1; j >= 0; j --) { sum = 0; for (k = 0; k < n; k ++) { if (k != j) { sum += buf1[j][k]; } } t = (single[j][c] - sum )/ buf1[j][j]; res[j][c] = t; for (i = j; i >= 1; i --) { buf1[i-1][j] *= t; } } return res; } float** reverse_matrix(float **buf, int n) { int i, j; float** buf_new = (float**)malloc(n * sizeof(float*)); for (i = 0; i < n; i ++) { buf_new[i] = (float*)malloc(n * sizeof(float)); } for (i = 0; i < n; i ++) { for (j = 0; j < n; j++) { buf_new[i][j] = buf[i][j]; } } float** single = (float**)malloc(n * sizeof(float*)); for (i = 0; i < n; i ++) { single[i] = (float*)malloc(n * sizeof(float)); } for (i = 0; i < n; i ++) { for (j = 0; j < n; j ++) { single[i][j] = 0; } single[i][i] = 1; } float** matrix; int k, t, p; p = 0; for(j = 0; j < n; j ++) { t = j + 1; k = 0; for(i = j; i < n - 1; i ++) { float x = buf[j][j]; float y = buf[t][j]; float cos = (x / sqrt(x*x + y*y)); float sin = - (y / sqrt(x*x + y*y)); matrix = matrix_return(cos, sin, n, j, k+p); buf = multiplication(matrix, buf, n); single = multiplication(matrix, single, n); t ++; k ++; } p ++; } float** result = (float**)malloc(n * sizeof(float*)); for (i = 0; i < n; i ++) { result[i] = (float*)malloc(n * sizeof(float)); } for(j = 0; j < n; j ++) { result = fill_matrix_for_reverse(buf, single, result, n, j); } return (result); } int main() { int i, j, N, K, k; float r; printf("%s\n", "Input N:"); scanf("%d", &N); printf("%s\n", "Input K:"); scanf("%d", &K); float t = 1.0/N; float h = 1.0/K; cout << "t: " << t << endl; cout << "h: " << h << endl; ofs << "N: " << N << endl; ofs << "K: " << K << endl; ofs << "t: " << t << endl; ofs << "h: " << h << endl; float *vector_u; vector_u = (float*)malloc(K*sizeof(float)); float *vector_v; vector_v = (float*)malloc(K*sizeof(float)); float *vector_b; vector_b = (float*)malloc(2*K*sizeof(float)); for (i = 1; i <= K; i ++) { vector_u[i] = (float)i/K; } vector_u[0] = 0; for (i = 0; i < K; i ++) { vector_v[i] = (float)(K-i)/K; } vector_v[K] = 0; //print_vector(vector_u, K+1); //print_vector(vector_v, K+1); float **matrix; matrix = (float**)malloc(2*K*sizeof(float*)); for(int i = 0; i < 2*K; i ++) { matrix[i] = (float*)malloc(2*K*sizeof(float)); } ofs << "Integral: " << endl; for (k = 0; k <= N; k ++) { float integral = 0; for (i = 1; i <= K-1; i ++) { integral += (vector_u[i]*vector_u[i]+vector_v[i]*vector_v[i]); } integral *= h; integral += (h/2)*(vector_u[K]*vector_u[K]+vector_v[0]*vector_v[0]); //ofs << "Iteration: " << k << "; Integral: " << integral << endl; ofs << integral << ", "; fill_matrix(matrix, 2*K, t, h, vector_u, vector_v, vector_b); //print_vector(vector_b, 2*K); //print_matrix(matrix, 2*K); float** rev_m = reverse_matrix(matrix, 2*K); //print_matrix(rev_m, 2*K); float* res_v = multiplication1(rev_m, vector_b, 2*K); //print_vector(res_v, 2*K); for (i = 1; i <= K; i ++) { vector_u[i] = res_v[i+K-1]; } vector_u[0] = 0; for (i = 0; i < K; i ++) { vector_v[i] = res_v[i]; } vector_v[K] = 0; } ofs << endl; ofs << "vector_u:" << endl; print_vector(vector_u, K+1); ofs << "vector_v:" << endl; print_vector(vector_v, K+1); return 0; }