lab6.cpp

The decision of a boundary problem for the equation of Puassona Method of splittings usual and повы

#include "StdAfx.h"
#include <iostream.h>
#include <conio.h>
#include <math.h>

const double a1 = 0;
const double a2 = 1;
const double b1 = 0;
const double b2 = 1;
const double eps = 0.0001;

#define nX 5
#define nY 5
#define hX (a2 - a1)/nX
#define hY (b2 - b1)/nY
#define m1 hX*hX/12
#define m2 hY*hY/12

const double t = 0.5*(hY*hX*hY*hY)/(hX*hX + hY*hY);

int i, j;
double u[nX+1][nY+1], v[nX+1][nY+1], z[nX][nY];

double Function_F(double x, double y)
{
	return 0;
}

void main()
{
	double norm, s1, s2, s3, count=0;
	
	for(i = 0; i <= nX; i++)
		for(j = 0; j <= nY; j++)
			v[i][j] = 0;
	
	for(i = 1; i < nY; i++)
	{
		u[0][i] = 0;
		u[nX][i] = 0;
	}
	for(i = 1; i < nX; i++)
	{
		u[i][nY] = 10000;
	}
	for(i = 0; i <= nX; i++)
	{
		u[i][0] = 0;
	}
	u[nX][nY] = 5000;
	u[0][nY] = 5000;
	
	for(i = 1; i < nX; i++)
		for(j = 1; j < nY; j++)
			u[i][j] = 0;

	do
	{
		s2 = 0;

		for(i = 1; i < nX; i++)
		{
			for(j = 1; j < nY; j++)
			{	
				z[i][j] = (u[i-1][j] - 2*u[i][j] + u[i+1][j])/(hX*hX) + 
					      (u[i][j-1] - 2*u[i][j] + u[i][j+1])/(hY*hY) -
						  (m1+m2+t)*(u[i-1][j-1] - 2*u[i][j-1] + u[i+1][j-1] - 
									 2*u[i-1][j] + 4*u[i][j] - 2*u[i+1][j] + 
									 u[i-1][j+1] - 2*u[i][j+1] + u[i+1][j+1])/(hX*hX*hY*hY) +
						  t*m1*(u[i-2][j-1] - 4*u[i-1][j-1] + 6*u[i][j-1] - 4*u[i+1][j-1] + u[i+2][j-1] - 
						     2*(u[i-2][j] - 4*u[i-1][j] + 6*u[i][j] - 4*u[i+1][j] + u[i+2][j]) + 
							    u[i-2][j+1] - 4*u[i-1][j+1] + 6*u[i][j+1] - 4*u[i+1][j+1] + u[i+2][j+1])/(hX*hX*hX*hX*hY*hY) + 
						  t*m2*(u[i-1][j-2] - 4*u[i-1][j-1] + 6*u[i-1][j] - 4*u[i-1][j+1] + u[i-1][j+2] -
						     2*(u[i][j-2] - 4*u[i][j-1] + 6*u[i][j] - 4*u[i][j+1] + u[i][j+2]) + 
							    u[i+1][j-2] - 4*u[i+1][j-1] + 6*u[i+1][j] - 4*u[i+1][j+1] + u[i+1][j+2])/(hX*hX*hY*hY*hY*hY) - 
						  t*m1*m2*(u[i-2][j-2] - 4*u[i-1][j-2] + 6*u[i][j-2] - 4*u[i+1][j-2] + u[i+2][j-2] - 
						        4*(u[i-2][j-1] - 4*u[i-1][j-1] + 6*u[i][j-1] - 4*u[i+1][j-1] + u[i+2][j-1]) + 
								6*(u[i-2][j] - 4*u[i-1][j] + 6*u[i][j] - 4*u[i+1][j] + u[i+2][j]) - 
								4*(u[i-2][j+1] - 4*u[i-1][j+1] + 6*u[i][j+1] - 4*u[i+1][j+1] + u[i+2][j+1]) + 
								   u[i-2][j+2] - 4*u[i-1][j+2] + 6*u[i][j+2] - 4*u[i+1][j+2] + u[i+2][j+2])/(hX*hX*hX*hX*hY*hY*hY*hY);
				
				u[i][j] += t*z[i][j];

				s1 = fabs((u[i][j] - v[i][j])*(u[i][j] - v[i][j]));
				s2 += s1;
			}
		}
		norm = sqrt(s2*hX*hY);

		for(i = 1; i < nX; i++)
			for(j = 1; j < nY; j++)
				v[i][j] = u[i][j];

		s3 = 0;
		for(i = 1; i < nX; i++)
			for(j = 1; j < nY; j++)
				if(norm >= eps)
					s3 = norm;
		count++;
	}while(s3 >= eps);

	cout<<"y/x";
	for(i = 0; i <= nX; i++)
		cout<<i*hX<<"    ";
	cout<<endl;
	for(j = 0; j <= nY; j++)
	{
		cout<<j*hY<<" ";
		for(i = 0; i <= nX; i++)
			cout<<u[i][j]<<" ";
		cout<<endl;
	}
}