microbik.ru
1 2 ... 5 6
Министерство образования и науки Российской Федерации

Новосибирский государственный технический университет

Кафедра прикладной математики


Лабораторная работа №1

по дисциплине «Уравнения математической физики»

РЕШЕНИЕ ЭЛЛИПТИЧЕСКИХ КРАЕВЫХ ЗАДАЧ

МЕТОДОМ КОНЕЧНЫХ РАЗНОСТЕЙ

Факультет: ПМИ

Группа: ПМ-86

Студент: Пустосмехов В.А.

Преподаватели: Персова М.Г.

Задорожный А.Г.

Токарева М.Г.

Вариант: 5

Новосибирск, 2011 г.

Цель работы

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

Область имеет L-образную область. Предусмотреть учет первых и третьих краевых условий.
Анализ задачи


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




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

Все узлы можно разделить на два вида. Краевые узлы – узлы, которые находятся на границе расчётной области. Регулярные узлы – узлы, которые находятся внутри расчётной области. Они различаются количеством соседних узлов и способом формирования строки матрицы.

По сути, метод заключается в том, что мы пробегаем по всем узлам сетки в порядке нумерации узлов и записываем соответствующие уравнения. Это может быть дискретный оператор Лапласа, если узел регулярный. Если узел принадлежит краевой границе, то в зависимости от типа краевого условия. Вместо исходной задачи получаем пятидиагональную СЛАУ, решение которой является значение искомой функции в узлах сетки.

Решение СЛАУ будем осуществлять с помощью метода Гаусса-Зейделя из лабораторной работы №2 по дисциплине «Численные методы».

Текст программы

Структура данных

struct LArea

{

double x1, y1, x2, y2, x3, x3;

double k1, k2, k3, k4;

int n1, n2, n3, n4;

double h01, h02, h03, h04;

};

Генератор сетки

#include

#include

#include

#include "LArea.h"

LArea A;

double** points; int pointsN;

double* stepsX; int M;

double* stepsY; int N;
bool input()

{

printf("----------- (x3,y3)\n");

printf("| |\n ");

printf("| |\n ");

printf("| | k3\n ");

printf("| |\n ");

printf("| |\n ");

printf("| | (x2,y2) k2\n");

printf("| k1 -------------------|\n");

printf("| |\n");

printf("| k4 |\n");

printf("| |\n");

printf("-----------------------------|\n");

printf("(x1,y1)\n\n");
printf("x1 = ");

scanf("%lf",&A.x1);
printf("y1 = ");

scanf("%lf",&A.y1);
printf("x2 = ");

scanf("%lf",&A.x2);
printf("y2 = ");

scanf("%lf",&A.y2);
printf("x3 = ");

scanf("%lf",&A.x3);
printf("y3 = ");

scanf("%lf",&A.y3);
printf("k1 = ");

scanf("%lf",&A.k1);
printf("k2 = ");

scanf("%lf",&A.k2);
printf("k3 = ");

scanf("%lf",&A.k3);
printf("k4 = ");

scanf("%lf",&A.k4);
printf("n1 = ");

scanf("%d",&A.n1);
printf("n2 = ");

scanf("%d",&A.n2);
printf("n3 = ");

scanf("%d",&A.n3);
printf("n4 = ");

scanf("%d",&A.n4);
if (A.x1 < A.x2 &&

A.x2 < A.x3 &&

A.y1 < A.y2 &&

A.y2 < A.y3 &&

A.k1 > 0 &&

A.k2 > 0 &&

A.k3 > 0 &&

A.k4 > 0 &&

A.n1 > 0 &&

A.n2 > 0 &&

A.n3 > 0 &&

A.n4 > 0)

return true;

else

return false;

}
void generateGrid()

{

int i, j;

double x, y, hx, hy, W;
pointsN = (A.n1 + 1) * (A.n3 + 1) // |

+ (A.n1 + A.n2 + 1) * (A.n4 + 1) - // _

- (A.n1 + 1); // На стыке узлы были посчитаны два раза
points = new double*[pointsN];

for (i=0; i < KolTop; i++)

Top[i] = new double[2];
if (A.k1 != 1.0)

A.h01 = (A.x2 - A.x1)*(1 - A.k1)/(1-pow(A.k1,A.n1));

else

A.h01 = (A.x2 - A.x1)/(A.n1);

if (A.k2 != 1.0)

A.h02 = (A.x3 - A.x2)*(1 - A.k2)/(1-pow(A.k2,A.n2));

else

A.h02 = (A.x3 - A.x2)/(A.n2);
if (A.k3 != 1.0)

A.h03 = (A.y3 - A.y2)*(1 - A.k3)/(1-pow(A.k3,A.n3));

else

A.h03 = (A.y3 - A.y2)/(A.n3);

if (A.k4 != 1.0)

A.h04 = (A.y2 - A.y1)*(1 - A.k4)/(1-pow(A.k4,A.n4));

else

A.h04 = (A.y2 - A.y1)/(A.n4);

// L

i = 0;

y = A.y1; hy = A.h04; W = hy / 10.0;

while (y <= A.y2 + W / 10000.0)

{

x = A.x1;

W = hx = A.h01;

while (x <= A.x2 + W/10000.0)

{

points[i][0] = x;

points[i][1] = y;

x += hx;

W = hx *= A.k1;

i++;

}

i += A.n2;

y += hy;

W = hy *= A.k4;

}
// _

i = A.n1 + 1;

hy = A.h04; y = A.y1; W = hy/10.0;

while (y <= A.y2 + W/10000.0)

{

hx = A.h02;

x = A.x2 + hx;

W = hx *= A.k2;

while (x <= A.x3 + W/10000.0)

{

points[i][0] = x;

points[i][1] = y;

x += hx;

W = hx *= A.k2;

i++;

}

y += hy;

W = hy *= A.k4;

i += A.n1 + 1;

}
// |

i = (A.n4 + 1)*(A.n1 + A.n2 +1);

hy = A.h03; y = A.y2 + hy;

hy *= A.k3; W = hy/10.0;

while (y <= A.y3 + W/10000.0)

{

x = A.x1;

W = hx = A.h01;

while (x <= A.x2 + W/10000.0)

{

points[i][0] = x;

points[i][1] = y;

x += hx;

W = hx *= A.k1;

i++;

}

y += hy;

W = hy *= A.k3;

}

}
void generateSteps()

{

int i;

double h, y, t;

M = A.n1 + A.n2;

N = A.n3 + A.n4;

stepsX = new double[M];

stepsY = new double[N];

for (i=0; i
stepsX[i] = points[i+1][0] - points[i][0];
y = A.y1;

h = A.h04;

t = A.y1;

y += h;

for(i=0; i
{

stepsY[i] = y - t;

t = y;

h *= A.k4;

y += h;

}
h = A.h03;

y = A.y2 + h;

t = A.y2;

for(; i
{

stepsY[i] = y - t;

t = y;

h *= A.k3;

y += h;

}

}
void output()

{

int i;

FILE* out = fopen("grid.txt", "w");
fprintf(out, "%d\n", pointsN);

for(i=0; i
fprintf(out,"%.15lf\t\t\t\t\t%.15lf\n",points[i][0],points[i][1]);

fclose(out);
out = fopen("steps.txt","w");

fprintf(out,"%d\n",M);

for(i=0; i
fprintf(out,"%.15lf\n",stepsX[i]);
fprintf(out,"\n\n%d\n",N);

for(i=0; i
fprintf(out,"%.15lf\n",stepsY[i]);

fclose(out);
out = fopen("n.txt","w");

fprintf(out,"%d\n",A.n1);

fprintf(out,"%d\n",A.n2);

fprintf(out,"%d\n",A.n3);

fprintf(out,"%d\n",A.n4);

fclose(out);

}
int main()

{

if (input())

{

generateGrid();

generateSteps();

output();

printf("Done!\n");

_getch();

return 0;

}

else

{

printf("?");

_getch();

return 1;

}

}

Генератор СЛАУ

#include

#include

#include

#include "LArea.h"

LArea A; // n1, n2, n3, n4

double** points; int pointsN;

double** A;

double* F;

int diOffsets[7]; // Смещения диагоналей

int borderTypes[6]; // Сторона -> Номер краевого условия

int** borderPoints; // Номера вершил на границе

double beta[6];
bool input()

{

int i, t;

FILE *in = fopen("grid.txt","r");

if (!in) return false;
fscanf(in,"%d",&pointsN);

points = new double*[pointsN];

F = new double[pointsN];

A = new double*[pointsN];

for(i=0; i

{

points[i] = new double[2];

A[i] = new double[7];

}

for(i=0; i

fscanf(in,"%lf %lf",&points[i][0],&points[i][1]);

fclose(in);
in = fopen("n.txt","r");

if (!in) return false;

fscanf(in,"%d",&A.n1);

fscanf(in,"%d",&A.n2);

fscanf(in,"%d",&A.n3);

fscanf(in,"%d",&A.n4);

fclose(in);
in = fopen("borderTypes.txt","r");

if(!in) return false;

for(i=0; i<6; i++){

fscanf(in,"%d",&t);

if(t==1 || t==3)

borderTypes[i] = t;

else return false;

}

fclose(in);
in = fopen("beta.txt","r");

if (!in) return false;

for(i=0; i<6; i++)

fscanf(in,"%lf",&beta[i]);

fclose(in);
return true;

}
double U(double x, double y)

{

return sin(x)+y;

}
double Ubeta(int t, double x, double y)

{

switch(t){

case 0: return (- + beta[0]*U(x,y))/beta[0]; break;

case 1: return (cos(x) + beta[1]*U(x,y))/beta[1]; break;

case 2: return (1 + beta[2]*U(x,y))/beta[2]; break;

case 3: return (cos(x) + beta[3]*U(x,y))/beta[3]; break;

case 4: return (1 + beta[4]*U(x,y))/beta[4]; break;

case 5: return (-cos(x) + beta[5]*U(x,y))/beta[5];break;

default: break;

}

}
double gamma(double x, double y)

{

return 1.0;

}
double f(double x, double y)

{

return sin(x) + gamma(x,y)*U(x,y);

}
void NullifyRow(int i)

{

for(int j=0; j<7; j++)

A[i][j] = 0.0;

}
// Номер границы если на границе, -1 если не на границе

int pointIsOnBorder(int t)

{

int i,j,n;

for(i=0; i<6; i++){

switch(i){

case 0: n = A.n1 + A.n2 + 1; break;

case 1: n = A.n4 + 1; break;

case 2: n = A.n2 + 1; break;

case 3: n = A.n3 + 1; break;

case 4: n = A.n1 + 1; break;

case 5: n = A.n3 + A.n4 + 1; break;

}

for(j=0; j
if(borderPoints[i][j] == t)

return i;

}

}

return (-1);

}
void fillBorderPoints()

{

int i,n,t;

//выделили память под 6 граней

borderPoints = new int*[6];
//выделили память под каждую грань в зависимости от кол-ва узлов в грани

for(i=0; i<6; i++){

switch(i){

case 0: n = A.n1 + A.n2 + 1; break;

case 1: n = A.n4 + 1; break;

case 2: n = A.n2 + 1; break;

case 3: n = A.n3 + 1; break;

case 4: n = A.n1 + 1; break;

case 5: n = A.n3 + A.n4 + 1; break;

}

borderPoints[i] = new int[n];

}
//заполняем вершинами

t = A.n1 + A.n2 + 1; //первая грань

for(i=0; i
borderPoints[0][i] = i;
t = A.n4 + 1; //вторая грань

borderPoints[1][0] = A.n1 + A.n2;

for(i=1; i
borderPoints[1][i] = borderPoints[1][i-1] + (A.n1 + A.n2 + 1);
t = A.n2 + 1; //третья грань

borderPoints[2][0] = borderPoints[1][A.n4];

for(i=1; i
borderPoints[2][i] = borderPoints[2][i-1]-1;
t = A.n3 +1; //четвертая грань

borderPoints[3][0] = borderPoints[2][A.n2];

borderPoints[3][1] = borderPoints[3][0] + A.n1 + A.n2 + 1;

for(i=2; i
borderPoints[3][i] = borderPoints[3][i-1] + A.n1 + 1;
t = A.n1 + 1; //пятая грань

borderPoints[4][0] = borderPoints[3][A.n3];

for(i=1; i
borderPoints[4][i] = borderPoints[4][i-1] - 1;
t = A.n3; //шестая грань, первая часть

borderPoints[5][0] = borderPoints[4][A.n1];

for(i=1; i
borderPoints[5][i] = borderPoints[5][i-1] - (A.n1 + 1);
i = t; t = A.n3 + A.n4 + 1; //шестая грань, вторая часть

for(; i
borderPoints[5][i] = borderPoints[5][i-1] - (A.n1 + A.n2 + 1);
//сформировали массив диагоналей I

diOffsets[3] = 0;
diOffsets[2] = -1;

diOffsets[4] = 1;
diOffsets[1] = -(A.n1 + 1);

diOffsets[5] = (A.n1 + 1);
diOffsets[0] = -(A.n1 + A.n2 + 1);

diOffsets[6] = (A.n1 + A.n2 + 1);

}
void fillRegularPoints(int i)

{

int t;

double h,h1;
//вычислили hi и hi-1 по Х

h = points[i+1][0] - points[i][0];

h1 = points[i][0] - points[i-1][0];
A[i][3] = -1.0*(2.0/(h*h1) + gamma(points[i][0],points[i][1]));//вычисляем U[l].первая часть (вторая часть по У считается)

A[i][2] = 2.0/(h1*(h+h1));//вычисляем U[l-1]

A[i][4] = 2.0/(h*(h+h1));//вычисляем U[l+1]
//вычислили hi и hi-1 по Y

t = (A.n1 + A.n2 + 1)*(A.n4 + 1);

if(i < t){//если узел в большОм нижнем прямоугольнике

h = points[i + A.n1+A.n2+1][1] - points[i][1];

h1 = points[i][1] - points[i - (A.n1+A.n2+1)][1];

A[i][0] = 2.0/(h1*(h+h1));//вычисляем U[l-kx]

A[i][6] = 2.0/(h*(h+h1));//вычисляем U[l+kx]

}

else{//иначе узел в верхнем прямоугольнике

if(i>=t&&(i<(t+A.n1))){

h = points[i + A.n1+1][1] - points[i][1];

h1 = points[i][1] - points[i - (A.n1+A.n2+1)][1];

A[i][0] = 2.0/(h1*(h+h1));//вычисляем U[l-kx]

A[i][5] = 2.0/(h*(h+h1));//вычисляем U[l+kx]

}

else{

h = points[i + A.n1+1][1] - points[i][1];

h1 = points[i][1] - points[i - (A.n1+1)][1];

A[i][1] = 2.0/(h1*(h+h1));//вычисляем U[l-kx]

A[i][5] = 2.0/(h*(h+h1));//вычисляем U[l+kx]

}

}
A[i][3] += -2.0/(h*h1);//вычисляем U[l].вторая часть

}
void FirstborderTypes(int i)

{

A[i][3] = 1.0;

F[i] = U(points[i][0],points[i][1]);

}
void Cornerpoints(int i, int p)

{

double h;

if(borderTypes[p+1] == 3 || !i){//если "смежная" граница тоже с 3м краевым, то производим сумму 3х краевых обоих границ

switch(p){

case 0:

if(i == 0){//на 0-й границе может быть 2 узла (0 и n1+n2+1)

if(borderTypes[0]==1 || borderTypes[5]==1)

FirstborderTypes(i);

else{//значит на обеих границах 3и краевые

//от 0й границы

h = points[i + (A.n1+A.n2+1)][1] - points[i][1];

A[i][3] = 1.0/h + beta[0];

A[i][6] = -1.0/h;

F[i] = beta[0]*Ubeta(0,points[i][0],points[i][1]);
//от 5й границы

h = points[i+1][0] - points[i][0];

A[i][3] += 1.0/h + beta[5];////!!!!!!!!!!!!!!!!!!!

A[i][4] += -1.0/h;

F[i] += beta[5]*Ubeta(5,points[i][0],points[i][1]);

}

}

else{//иначе это не 0й узел, а n1+n2+1-й

//от 0й границы

h = points[i + (A.n1+A.n2+1)][1] - points[i][1];

A[i][3] = 1.0/h + beta[0];

A[i][6] = -1.0/h;

F[i] = beta[0]*Ubeta(0,points[i][0],points[i][1]);
//от 1границы

h = points[i][0] - points[i-1][0];

A[i][3] += 1/h + beta[1];

A[i][2] += -1.0/h;

F[i] += beta[1]*Ubeta(1,points[i][0],points[i][1]);

}

break;

case 1:

//от 1границы

h = points[i][0] - points[i-1][0];

A[i][3] = 1/h + beta[1];

A[i][2] = -1.0/h;

F[i] = beta[1]*Ubeta(1,points[i][0],points[i][1]);
//от 2й границы

h = points[i][1] - points[i-(A.n1+A.n2+1)][1];

A[i][3] += 1.0/h + beta[2];

A[i][0] += -1.0/h;

F[i] += beta[2]*Ubeta(2,points[i][0],points[i][1]);

break;

case 2:

//от 2й границы

h = points[i][1] - points[i-(A.n1+A.n2+1)][1];

A[i][3] = 1.0/h + beta[2];

A[i][0] = -1.0/h;

F[i] = beta[2]*Ubeta(2,points[i][0],points[i][1]);
//от 3й границы

h = points[i][0] - points[i-1][0];

A[i][3] += 1/h + beta[3];

A[i][2] += -1.0/h;

F[i] += beta[3]*Ubeta(3,points[i][0],points[i][1]);

break;

case 3:

//от 3й границы

h = points[i][0] - points[i-1][0];

A[i][3] = 1/h + beta[3];

A[i][2] = -1.0/h;

F[i] = beta[3]*Ubeta(3,points[i][0],points[i][1]);
//от 4й границы

h = points[i][1] - points[i-(A.n1+1)][1];

A[i][3] += 1/h + beta[4];

A[i][1] += -1.0/h;

F[i] += beta[4]*Ubeta(4,points[i][0],points[i][1]);

break;

case 4:

//от 4й границы

h = points[i][1] - points[i-(A.n1+1)][1];

A[i][3] = 1/h + beta[4];

A[i][1] = -1.0/h;

F[i] = beta[4]*Ubeta(4,points[i][0],points[i][1]);
//от 5й границы

h = points[i+1][0] - points[i][0];

A[i][3] += 1.0/h + beta[5];///!!!!!!!!!!!

A[i][4] += -1.0/h;

F[i] += beta[5]*Ubeta(5,points[i][0],points[i][1]);

break;

default:

break;

}

}

else//иначе значит в соседней области 1е краевые. вычисляем их,они приоритетнее

FirstborderTypes(i);

}
void ThirdborderTypes(int i, int p)

{

int t1,t2,t3,t4,t5;

double h;

t1 = A.n1 + A.n2; //номер 2го углового узла

t2 = (t1+1)*(A.n4 + 1)-1;//номер 3го углового узла

t3 = t2 - A.n2; //номер 4го углового узла

t4 = pointsN - 1; //номер 5го углового узла

t5 = t4 - A.n1; //номер 6го углового узла

if(i==0 || i==t1 || i==t2 || i==t3 || i==t4 || i==t5)//если краевая точка угловая (принадлежит 2м вершинам)

Cornerpoints(i,p);

else{//если краевая точка строго принадлежит одной грани

switch(p){

case 0:

h = points[i + (A.n1+A.n2+1)][1] - points[i][1];
A[i][3] = 1.0/h + beta[0];

A[i][6] = -1.0/h;

F[i] = beta[0]*Ubeta(0,points[i][0],points[i][1]);

break;

case 1:

h = points[i][0] - points[i-1][0];
A[i][3] = 1/h + beta[1];

A[i][2] = -1.0/h;

F[i] = beta[1]*Ubeta(1,points[i][0],points[i][1]);

break;

case 2:

h = points[i][1] - points[i-(A.n1+A.n2+1)][1];
A[i][3] = 1.0/h + beta[2];

A[i][0] = -1.0/h;

F[i] = beta[2]*Ubeta(2,points[i][0],points[i][1]);

break;

case 3:

h = points[i][0] - points[i-1][0];
A[i][3] = 1/h + beta[3];

A[i][2] = -1.0/h;

F[i] = beta[3]*Ubeta(3,points[i][0],points[i][1]);

break;

case 4:

h = points[i][1] - points[i-(A.n1+1)][1];
A[i][3] = 1/h + beta[4];

A[i][1] = -1.0/h;

F[i] = beta[4]*Ubeta(4,points[i][0],points[i][1]);

break;

case 5:

h = points[i+1][0] - points[i][0];
A[i][3] = 1.0/h + beta[5];

A[i][4] = -1.0/h;

F[i] = beta[5]*Ubeta(5,points[i][0],points[i][1]);

break;

default:

break;

}

}

}
void EDGECondition(int i, int p)//обработка вершины,если это краевоя вершина

{

switch(borderTypes[p]){

case 1: FirstborderTypes(i); break;

//case 2: SecondborderTypes(i); break;

case 3: ThirdborderTypes(i,p); break;

default: break;

}

}
void MethodFinalDifferences()

{

int i,p;

CreateborderPoints(); //создаем массив с вершинами,которые на границе

for(i=0; i

printf("\t%d / %d\n",i+1,pointsN);

NullifyRow(i); //обнуляем строку в матрице

F[i] = -f(points[i][0],points[i][1]); //определили значение вектора правой части

if((p=pointIsOnBorder(i))==-1) //если i-я вершина не на границе

fillRegularPoints(i); //то обрабатываем ее как регулярную

else

EDGECondition(i,p); //иначе - как краевую

}

}
void Output()

{

int i,j;

FILE *out = fopen("output.txt","w");
fprintf(out,"%d\n1.e-21\n10000\n\n",pointsN);
for(i=0; i

for(j=0; j<7; j++){

fprintf(out,"%.15lf\t\t",A[i][j]);

}

fprintf(out,"\n");

}
fprintf(out,"\n");

for(i=0; i<7; i++)

fprintf(out,"%d ",diOffsets[i]);
fprintf(out,"\n\n");

for(i=0; i

fprintf(out,"%.15lf\n",F[i]);
fprintf(out,"\n\n");

for(i=0; i

fprintf(out,"0 ");
fclose(out);

}
void Work()

{

int i,j;

double sum;

for(i=0; i

sum = 0.0;

for(j=0; j<7; j++)

sum += A[i][j];

printf("\t%.15lf\n",sum);

}

}
int main()

{

if(input()){

MethodFinalDifferences();

Output();

//Work();

printf("\n\n\tComplete!");

_getch();

return 0;

}

else{

printf("\n\n\tError in input data!!!");

_getch();

return 1;

}

}


следующая страница >>