Задача Дирихле - (реферат)
Дата добавления: март 2006г.
Задача Дирихле
1. ПОСТАНОВКА ЗАДАЧИ
Решить численно задачу Дирихле для уравнения Лапласа :
(x, y)? ?D; u|Г=xy2=f(x, y) ;
область D ограничена линиями: x=2 , x=4 , y=x , y=x+4 ;
(x0, y0 ) = (3, 5) .
Следует решить задачу сначала с шагом по x и по y : h=0. 2, потом с шагом h=0. 1 . Точность решения СЛАУ? =0. 01 .
2. ОПИСАНИЕ МЕТОДА РЕШЕНИЯ ПОСТАВЛЕННОЙ ЗАДАЧИ
Поставленная задача решается численно с помощью программы, реализующей метод сеток , разработанный для численного решения задачи Дирихле для уравнений эллептического типа.
Программа написана на языке C++ , в среде Borland C++ версии 3. 1. Ниже описан алгоритм работы этой программы.
1. На первом шаге область D дискретизируется. Она заменяется на область Dh путем разбиения области D параллельными прямыми по следующему правилу: yi=y0 ? ih, xj=x0 ? ih , i, j=0, 1, 2…. PР. Разбиение производится до тех пор, пока текущая прямая не будет лежать целиком вне области D. Получается множество точек (xi, yj). 2. За область Dh принимают те точки множества (xi, yj) , которые попали внутрь области D, а также дополняют это множество граничными точками.
3. Во всех точках области Dh вычисляются значения функции f(xi, yj) . 4. За область Dh* принимаются все внутренние точки области Dh, т. е. удовлетворяющие требованию: (xi, yj)? ? Dh* , если (xi+1, yj) ? Dh , (xi-1, yj)? ? Dh , (xi, yj+1)? ? Dh , (xi, yj-1)? ? Dh . 5. Во всех точках области Dh* вычисляется функция F(N)*[i, j] ( индекс N обозначает номер итерации, на которой производится вычисление):
F(N)*[i, j]=(f(xi+1, yj) + f(xi-1, yj) + f(xi, yj+1)? ?? f(xi, yj-1))/4 6. Теперь если max | F(N+1)*[i, j] - F(N)*[i, j]|
3. ТЕКСТ ПРОГРАММЫ
#include
#include
#include
#include
#include
int i, j, k; // Variables
float h, x, y, tmp, E1;
struct point {
float xx;
float yy;
int BelongsToDh_;
int BelongsToDh;
float F;
float F_;
} p0, arrayP[13][33];
float arrayX[13];
float arrayY[33];
float diff[500];
void CreateNet(void); // Procedure Prototypes int IsLineFit(float Param);
void CrMtrD(void);
void RegArrayX();
void RegArrayY();
void CreateDh_();
int IsFit(point Par);
void FillF();
void CreateDh();
int IsInner(int i, int j);
void FillF_();
void CountDif();
void MakeFile();
void main(void) //MAIN
{
clrscr();
p0. xx = 3;
p0. yy = 5;
h = 0. 2;
p0. BelongsToDh_=1;
p0. BelongsToDh=1;
CreateNet();
RegArrayX();
RegArrayY();
CrMtrD();
CreateDh_();
FillF();
CreateDh();
FillF_();
CountDif();
while (E1>=0. 005) {
for(i=0; i for(j=0; j FillF_();
CountDif();
}
cout MakeFile();
getchar();
} //MAIN END
int IsLineFit(float par, char Axis) // does the line belong to the defined area {
switch(Axis) (par else return 0;
case 'x': if (par else if (par>4. 0) return 1;
else return 0;
}
void CreateNet(void) // Creation of Net (area D )
{
x = p0. xx;
i=0;
arrayX[i]=x;
while (! IsLineFit(x, 'x'))
{
x += h;
i++;
arrayX[i] = x;
}
x = p0. xx-h;
i++;
arrayX[i]=x;
while (! IsLineFit(x, 'x'))
{
x -= h;
i++;
arrayX[i] = x;
}
for (i=0; i printf("\n");
y = p0. yy;
i = 0;
arrayY[i]=y;
while (! IsLineFit(y, 'y'))
{
y += h;
i++;
arrayY[i] = y;
}
y = p0. yy - h;
i++;
arrayY[i]=y;
while (! IsLineFit(y, 'y'))
{
y -= h;
i++;
arrayY[i] = y;
}
for(i=0; i printf("\n");
} // end CreateNet
void RegArrayX() // Regulation of arrays X & Y
{
int LastUnreg = 13 ;
while (LastUnreg ! = 0) {
for(i=0; i if (arrayX[i]>arrayX[i+1]) {double tmp=arrayX[i];
arrayX[i]=arrayX[i+1];
arrayX[i+1]=tmp; }}
LastUnreg=LastUnreg-1; }
for (i=0; i } printf("\n");
}
void RegArrayY()
{
int LastUnreg = 33 ;
while (LastUnreg ! = 0) {
for(i=0; i if (arrayY[i]>arrayY[i+1]) { tmp=arrayY[i];
arrayY[i]=arrayY[i+1];
arrayY[i+1]=tmp; }}
LastUnreg=LastUnreg-1; }
for (i=0; i printf("\n"); } // End of Regulation
void CrMtrD(void) //Create general Matrix
{
for(i=0; i for(j=0; j arrayP[i][j]. BelongsToDh=0; }
for(i=0; i for(j=0; j arrayP[i][j]. xx=arrayX[i];
arrayP[i][j]. yy=arrayY[j];
}
// printf("%g %g", arrayP[12][0]. xx, arrayP[12][0]. yy);
// printf("\n");
}
int IsFit(point Par) //does point belong to area D?
{
if ((Par. xx=1. 99) && (Par. yy>=Par. xx)
&& (Par. yy else return 0;
}
void CreateDh_(void) //Create area Dh_
{
for(i=0; i for(j=0; j if (IsFit(arrayP[i][j])) arrayP[i][j]. BelongsToDh_=1;
cout cout }
void FillF(void) // calc function F(x, y) at area Dh_
{
for(i=0; i for(j=0; j if (arrayP[i][j]. BelongsToDh_==1)
arrayP[i][j]. F=arrayP[i][j]. xx*pow(arrayP[i][j]. yy, 2);
else arrayP[i][j]. F=0;
}
int IsInner(int i, int j) //Is point inner?
{
if ((arrayP[i-1][j]. BelongsToDh_==1) &&
(arrayP[i+1][j]. BelongsToDh_==1) &&
(arrayP[i][j+1]. BelongsToDh_==1) &&
(arrayP[i][j-1]. BelongsToDh_==1)) return 1;
else return 0;
}
void CreateDh(void) //Create area Dh
{
for(i=0; i for(j=0; j if ((arrayP[i][j]. BelongsToDh_==1) &&
IsInner(i, j))
arrayP[i][j]. BelongsToDh=1;
}
void FillF_() //calc new appr. values of F
{
for(i=0; i for(j=0; j if (arrayP[i][j]. BelongsToDh==1)
arrayP[i][j]. F_=(arrayP[i-1][j]. F+arrayP[i+1][j]. F+
arrayP[i][j-1]. F+arrayP[i][j+1]. F)/4;
else arrayP[i][j]. F_=0; }
}
void CountDif() // find maximal difference abs(F-F_)
{
k=0;
for(i=0; i for(j=0; j { if (arrayP[i][j]. BelongsToDh==1) {
diff[k]=fabs(arrayP[i][j]. F_-arrayP[i][j]. F);
k++; }}
E1=diff[0];
for (k=1; k if (diff[k]>E1) E1=diff[k]; }
}
void MakeFile()
{
ofstream f;
FILE *f1=fopen("surf. dat", "w1");
fclose(f1);
f. open("surf. dat", ios: :out, 0);
for(i=0; i for(j=0; j f " " f. close() ;
}
4. ГРАФИКИ РЕШЕНИЙ
РИС. 1 шаг h=0. 2
РИС. 2 шаг h=0. 1
5. ВЫВОД
Функция f(x, y) является неотрицательной в области D. Полученное решение лежит целиком над плоскостью XOY . Для данного решения выполняется принцип максимума.