.

Решение задачи Дирихле для уравнения Лапласа методом сеток

Язык: русский
Формат: реферат
Тип документа: Word Doc
0 659
Скачать документ

ПОСТАНОВКА ЗАДАЧИ

Решить численно задачу Дирихле для уравнения Лапласа :

(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 .

ОПИСАНИЕ МЕТОДА РЕШЕНИЯ ПОСТАВЛЕННОЙ ЗАДАЧИ

Поставленная задача решается численно с помощью программы,
реализующей метод сеток , разработанный для численного решения задачи
Дирихле для уравнений эллептического типа.

Программа написана на языке C++ , в среде Borland C++ версии 3.1. Ниже
описан алгоритм работы этой программы.

1. На первом шаге область D дискретизируется. Она заменяется на область
Dh путем разбиения области D параллельными прямыми по следующему
правилу: yi=y0 ? ih, xj=x0 ? ih , i,j=0,1,2…. РР
Разбиение производится до тех пор, пока текущая прямая не будет лежать
целиком вне области 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]|

#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<13;i++)for(j=0;j<33;j++) arrayP[i][j].F=arrayP[i][j].F_;FillF_();CountDif();}cout<<(0-arrayP[7][17].F_);MakeFile();getchar();} //MAIN ENDint IsLineFit(float par,char Axis) // does the line belong to the defined area{switch(Axis) {case 'y': if ((par>8.0) || (par<2.0)) return 1;else return 0;case 'x': if (par<1.9) return 1;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<13;i++) { printf("%g ",arrayX[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<33;i++) { printf("%g ",arrayY[i]);}printf("\n");} // end CreateNetvoid RegArrayX() // Regulation of arrays X & Y{int LastUnreg = 13 ;while (LastUnreg != 0) {for(i=0;iarrayX[i+1]) {double tmp=arrayX[i];

arrayX[i]=arrayX[i+1];

arrayX[i+1]=tmp;}}

LastUnreg=LastUnreg-1; }

for (i=0;i<13;i++) { printf("%g ",arrayX[i]);} printf("\n");}void RegArrayY(){int LastUnreg = 33 ;while (LastUnreg != 0) {for(i=0;iarrayY[i+1]) { tmp=arrayY[i];

arrayY[i]=arrayY[i+1];

arrayY[i+1]=tmp;}}

LastUnreg=LastUnreg-1; }

for (i=0;i<33;i++) { printf("%g ",arrayY[i]); }printf("\n");} // End of Regulationvoid CrMtrD(void) //Create general Matrix{for(i=0;i<13;i++)for(j=0;j<33;j++) {arrayP[i][j].BelongsToDh_=0;arrayP[i][j].BelongsToDh=0;}for(i=0;i<13;i++)for(j=0;j<33;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<=4) && (Par.xx>=1.99) && (Par.yy>=Par.xx)

&& (Par.yy<=Par.xx+4)) return 1;else return 0;}void CreateDh_(void) //Create area Dh_{for(i=0;i<13;i++)for(j=0;j<33;j++)if (IsFit(arrayP[i][j])) arrayP[i][j].BelongsToDh_=1;cout << arrayP[1][1].BelongsToDh_<< "\n";cout << arrayP[1][1].xx << " " << arrayP[1][1].yy<<"\n";}void FillF(void) // calc function F(x,y) at area Dh_{for(i=0;i<13;i++)for(j=0;j<33;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<13;i++)for(j=0;j<33;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<13;i++)for(j=0;j<33;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<13;i++)for(j=0;j<33;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<500;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<13;i++) for(j=0;j<33;j++?????????????????????????????????????????????•?????????? ?????????????†††?????ГРАФИКИ РЕШЕНИЙРИС.1 шаг h=0.2РИС.2 шаг h=0.15.ВЫВОДФункция f(x,y) является неотрицательной в области D. Полученное решение лежит целиком над плоскостью XOY . Для данного решения выполняется принцип максимума.PAGEPAGE 9

Нашли опечатку? Выделите и нажмите CTRL+Enter

Похожие документы
Обсуждение

Оставить комментарий

avatar
  Подписаться  
Уведомление о
Заказать реферат
UkrReferat.com. Всі права захищені. 2000-2019