Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <iomanip>
- #include <fstream>
- #include <stdlib.h>
- using namespace std;
- const int A = 3;
- const int N1 = 20;
- const int N2 = 20;
- const double phi = 0;
- const double Eps = 0.001;
- const double l1 = 1;
- const double l2 = 1;
- const double h1 = l1 / N1;
- const double h2 = l2 / N2;
- double k1(double x1, double x2) {
- return (x1 * x1 + x2 * x2 + 10);
- }
- double k2(double x1, double x2) {
- return (5*(x1 * x1 + x2 * x2) + 20);
- }
- double p1(double x1, double x2) {
- return 0/*(-2*(x1 + 5*x2))*/;
- }
- double p2(double x1, double x2) {
- return 0/*(-1*(x1 * x1 + 5*x2*x2))*/;
- }
- double q(double x1, double x2) {
- return 0 /*(5*x1*x1 + x2*x2 + 30)*/;
- }
- double func(double x1, double x2){
- return(A * exp(x1 + x2));
- }
- double Af(double *mas[], double *x[] , int i , int j) {
- double T = (k1(x[0][i] + h1/2, x[1][j]) * (mas[i+1][j] - mas[i][j])) / (h1*h1)-
- (k1(x[0][i] - h1/2, x[1][j]) * (mas[i][j] - mas[i-1][j])) / (h1*h1)+
- (k2(x[0][i], x[1][j] + h2/2) * (mas[i][j+1] - mas[i][j])) / (h2*h2)-
- (k2(x[0][i], x[1][j] - h2/2) * (mas[i][j] - mas[i][j-1])) / (h2*h2)+
- (p1(x[0][i], x[1][j]) * (mas[i+1][j] - mas[i-1][j])) / (2*h1)+
- (p2(x[0][i], x[1][j]) * (mas[i][j+1] - mas[i][j-1])) / (2*h2)-
- q(x[0][i], x[1][j]) * mas[i][j];
- cout << T<< endl;
- system("pause");
- return T;
- }
- int main() {
- setlocale(LC_ALL, "");
- double** x = new double*[2];
- x[0] = new double[N1 + 1];
- x[1] = new double[N2 + 1];
- for (int i(0); i <= N1; i++) {
- x[0][i] = i * h1;
- }
- for (int j(0); j <= N2; j++) {
- x[1][j] = j * h2;
- }
- double** y = new double*[N1+1];
- double** u = new double*[N1+1];
- for(int i(0); i<=N1; i++) {
- y[i] = new double[N2+1];
- u[i] = new double[N2+1];
- for(int j(0); j<=N2; j++) {
- y[i][j] = 0;
- u[i][j] = 0;
- }
- }
- for(int i(0); i<=N1; i++) {
- for(int j(1); j<N2; j++) {
- y[i][0] = A * exp(x[0][i]);
- y[i][N2] = A * exp(l2+x[0][i]);
- y[0][j] = A * exp(x[1][j]);
- y[N1][j] = A * exp(l1+x[1][j]);
- u[i][0] = A * exp(x[0][i]);
- u[i][N2] = A * exp(l2+x[0][i]);
- u[0][j] = A * exp(x[1][j]);
- u[N1][j] = A * exp(l1+x[1][j]);
- }
- }
- int count = 0;
- double norm;
- do {
- count++;
- double** r = new double* [N1];
- for(int i(0); i<=N1; i++) {
- r[i] = new double[N2];
- }
- for(int i(1); i<N1; i++) {
- for(int j(1); j<N2; j++) {
- r[i][j] = Af(u, x, i, j) - phi;
- }
- }
- double scalR = 0;
- for(int i(1); i<N1; i++) {
- for(int j(1); j<N2; j++) {
- scalR += r[i][j] * r[i][j] * h1*h1;
- }
- }
- double scalA = 0;
- for(int i(1); i<N1; i++) {
- for(int j(1); j<N2; j++) {
- scalA += Af(r, x, i, j) * r[i][j] * h1*h1;
- }
- }
- double tau = scalR / scalA;
- //Y[i][j]
- for(int i(1); i<N1; i++) {
- for(int j(1); j<N2; j++) {
- y[i][j] = u[i][j] - r[i][j] * tau;
- }
- }
- //Norm
- norm = -1;
- for(int i(1); i<N1; i++) {
- for(int j(1); j<N2; j++) {
- if(fabs(y[i][j]-u[i][j]) > norm) {
- norm = fabs(y[i][j]-u[i][j]);
- }
- }
- }
- //
- for(int i(1); i<N1; i++) {
- for(int j(1); j<N2; j++) {
- u[i][j]=y[i][j];
- }
- }
- } while(norm > Eps);
- for(int i(0); i<=N1; i++) {
- for(int j(0); j<=N2; j++) {
- cout<< setw(5) << setprecision(3) << func(x[0][i], x[1][j])<< " ";
- }
- cout<<endl<<endl;
- }
- cout << endl;
- cout<<"Count: " << count <<endl;
- for(int i(0); i<=N1; i++) {
- for(int j(0); j<=N2; j++) {
- cout<< setw(5) << setprecision(3) << y[i][j]<< " ";
- }
- cout<<endl<<endl;
- }
- system("pause");
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement