Решение задачи Неймана для уравнения Пуассона в прямоугольной области

Курсовая работа

«Решение задачи Неймана для уравнения Пуассона в прямоугольной области»

Содержание Введение Аннотация

1. Численная постановка задачи

2. Решение заданного примера

3. Листинг программы

4. Результаты работы программы Литература

Введение

Актуальность решения уравнения Пуассона состоит в том что, при вычислении решения системы уравнений Навье-Стокса, алгоритм решения этой системы происходит в несколько этапов одним из сложных является решение уравнения Пуассона.

Аннотация Разработать программу для решения дифференциального уравнения Лапласа:

в прямоугольной области ABCD, с вершинами A (0;0), B (0:1), C (1;1), D (1,0), шагом h=l/n, где n-количество узлов, принимающее условия Дирихле на всех границах кроме правой (на правой поставлено условие Неймана). Уравнение решается методом сеток, c точностью е= 0,1. Программа разработана на языке C++.

• 1. Численная постановка задачи

Уравнение Лапласа является модельным для эллиптических уравнений в частных производных.

Некоторые важные задачи, часто встречающиеся в приложениях, сводятся к решению одного эллиптического уравнения. К ним относятся задачи расчета дозвукового безвихревого (потенциального) течения газа и определения стационарного поля температуры в твердом теле.

В данной работе требуется решить, конечно — разностную задачу Дирихле или Неймана для уравнения Лапласа в прямоугольной области т. е. найти непрерывную функцию u (х, у), удовлетворяющую внутри многоугольной области

уравнению Лапласа

(1)

и принимающую на границе области заданные значения, т. е.

,

,

где f₁=0, f₂=0, f₃=0, f₄=0.

Будем считать, что u (х, у) непрерывна на границе области, т. е.

.

Выбрав шаги h, l по x и y соответственно, строим сетку:

, , ,

где, .

Вводя обозначения

аппроксимируем частные производные и в каждом внутреннем узле сетки центральными разностными производными второго порядка

и заменим уравнение Лапласа конечно-разностным уравнением

(2)

.

Погрешность замены дифференциального уравнения разностным, составляет величину .

Уравнения (2) вместе со значениями в граничных узлах образуют систему линейных алгебраических уравнений относительно приближенных значений функции и (х, у) в узлах сетки. Наиболее простой вид имеет эта система при :

, ,, (3)

.

При получении сеточных уравнений (3) была использована схема узлов. Набор узлов, используемых для аппроксимации уравнения в точке, называется шаблоном. В данной работе используется шаблон типа «крест».

Численное решение задачи Дирихле для уравнения Лапласа в многоугольнике состоит в нахождении приближенных значений искомой функции u (х, у) во внутренних узлах сетки. Для определения величин требуется решить систему линейных алгебраических уравнений (3).

В данной работе она решается методом минимальных невязок, который состоит в построении последовательности итераций вида:

где — невязка

— итерационный параметр

(верхним индексом s обозначен номер итерации). При последовательность сходится к точному решению системы (3). В качестве условия окончания итерационного процесса можно принять

.

Таким образом, погрешность приближенного решения, полученного методом сеток, складывается из двух погрешностей: погрешности аппроксимации дифференциального уравнения разностными значениями; погрешности, возникающей в результате приближенного решения системы разностных уравнений (3).

Известно, что описанная здесь разностная схема обладает свойством устойчивости и сходимости. Устойчивость схемы означает, что малые изменения в начальных данных приводят к малым изменениям решения разностной задачи. Только такие схемы имеет смысл применять в реальных вычислениях. Сходимость схемы означает, что при стремлении шага сетки к нулю () решение разностной задачи стремится к решению исходной задачи. Таким образом, выбрав достаточно малый шаг h, можно как угодно точно решить исходную задачу.

Пример решения такой задачи Дирихле, приведен в решении заданного примера.

2. Решение заданного примера программа уравнение лаплас прямоугольный Используя метод сеток, составить приближенное решение задачи Дирихле для уравнения Лапласа (1).

Решение получить в квадрате ABCD, с вершинами A (0;0), B (0:1), C (1;1), D (1,0) шагом h=l/n, где n-количество узлов, принимающее условия Дирихле на всех границах кроме правой (на правой границе поставлено условие Неймана).

;;; .

Систему линейных алгебраических уравнений решить по методом минимальных невязок, при е=0,1.

1) Построим сетку с шагом h=l=0,2

2. Построим итерационный процесс В виде начального приближения возьмем ,

условия окончания итерационного процесса: .

3. Листинг программы

#define _USE_MATH_DEFINES

#include

#include

#include

#include

#include

using namespace std;

//объявление типов точек области начало

#define FictivePoint 0 //не расчетная точка

#define ActualPoint 1 //расчетная точка

#define TopPoint 2 //точка на верхней границе

#define RightPoint 3 //точка на правой границе

#define BottomPoint 4 //точка на нижней границе

#define LeftPoint 5 //точка на левой границе

#define LeftBottomPoint 6 //точка на выпуклом угле

//объявление типов точек области конец

FILE *file;

//проверка на симетричность начало

double Simetric (int n, double **M)

{

int k=0;

for (int i=0;i

{

for (int j=0;j

{

if (M[i][j]==M[j][i])

{

k++;

}

}

}

if (k==n*n)

{

for (int i=0;i

{

for (int j=0;j

{

if (M[i][j]==M[j][i])

{

if (i≠j)

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 2) ;

}

if (i==j)

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

}

cout<<<" «;

}

}

cout<

}

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 2) ;

cout<<" Матрица симетрична!" <

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

}

if (k≠n*n)

{

for (int i=0;i

{

for (int j=0;j

{

if (M[i][j]==M[j][i])

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

cout<<<" «;

}

if (M[i][j]≠M[j][i])

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 4) ;

cout<<<" «;

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

}

}

cout<

}

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 4) ;

cout<<" Матрица не симетрична!" <

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

}

return 0;

}

//проверка на симетричность конец

//проверка на знакоопределенность начало

double More_then₀(int n, double **M)

{

double d;

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 9) ;

cout<<" Приводим матрицу к треугольному виду…" <

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

for (int k=0;k

{

for (int i=k+1;i

{

if (M[k][k]==0)

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 4) ;

printf («Произошло деление на 0! nДля выхода нажмите любую клавишу… n»);

_getch ();

exit (0);

}

else

{

d = M[i][k]/M[k][k];

for (int j=0;j

{

M[i][j]=M[i][j]-M[k][j]*d;

}

}

}

}

int z=0,p=0;

double DetM=1;

double *detM;

detM=new double[n];

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 9) ;

cout<<" Находим главные миноры…" <

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

for (int i=0;i

{

DetM*=M[i][i];

detM[i]=DetM;

printf («минор[%d] = %lfn», i+1,detM[i]);

}

for (int i=0;i

{

if (detM[i]>0)

{

z++;

}

}

if (z==n)

{

cout<<" Так как все главные миноры имеют положительный знак то: «<

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 2) ;

cout<<" Матрица положительно определена!" <

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

}

else

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 4) ;

cout<<" Матрица не знакоопределена!" <

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 7) ;

}

return 0;

}

//проверка на знакоопределенность конец

//оператор Пуассона начало

void A (int n, int **mask, double **u, double **res_A)

{

double h=1./n;

for (int i=0; i<=n; i++)

{

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

{

if (mask[i][j]==ActualPoint)

{

res_A[i][j]= - ((u[i+1][j]-2*u[i][j]+u[i-1][j]+u[i][j+1]-2*u[i][j]+u[i][j-1])/(h*h));

}

if (mask[i][j]==RightPoint)

{

res_A[i][j]= (u[i][j]-u[i-1][j])/(h*h);

}

if (mask[i][j]==TopPoint)

{

res_A[i][j]= (u[i][j]-u[i][j-1])/(h*h);

}

if (mask[i][j]==BottomPoint)

{

res_A[i][j]= (u[i][j]-u[i][j+1])/(h*h);

}

if (mask[i][j]==LeftPoint)

{

res_A[i][j]= (u[i][j]-u[i+1][j])/(h*h);

}

if (mask[i][j]==LeftBottomPoint)

{

res_A[i][j]=-(u[i+1][j]-4*u[i][j]+u[i][j+1])/(h*h); }

}

}

//оператор Пуассона конец

//разность векторов начало

void Substruction (int n, int **mask, double **a, double **b, double **res_S)

{

for (int i=0; i<=n; i++)

{

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

{

if (mask[i][j]≠FictivePoint)

{

res_S[i][j]=a[i][j]-b[i][j];

}

}

}

}

//разность векторов конец

//умножение вектора на скаляр начало

double Scalar (int n, int **mask, double **a, double **b)

{

double tmp=0;

for (int i=0;i<=n;i++)

{

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

{

if (mask[i][j]≠FictivePoint)

{

tmp=tmp+a[i][j]*b[i][j];

}

}

}

return tmp;

}

//умножение вектора на скаляр конец

//вычисление нормы вектора начало

double Norm (int n, int **mask, double **a)

{

double tmp=0;

for (int i=0;i<=n;i++)

{

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

{

if (mask[i][j]≠FictivePoint)

{

tmp=tmp+a[i][j]*a[i][j];

}

}

}

tmp=sqrt (tmp);

return tmp;

}

//вычисление нормы вектора конец

//вывод начало

void print (int n, double **a, double *f)

{

fprintf (file," n");

for (int i=0; i<=n-1; i++)

{

for (int j=0; j<=n-1; j++)

{

fprintf (file," %5.0lf «, a[i][j]);

}

fprintf (file," | %5.0lfn", f[i]);

}

fprintf (file," n");

}

//вывод конец

//копирование вектора в вектор начало

double copy (int n, double **a, double **b)

{

for (int i=0;i<=n;i++)

{

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

{

b[i][j]=a[i][j];

}

}

return 0;

}

//копирование вектора в вектор конец

int main ()

{

int n=0,**mask, count_it=0;

double **u,**u_accur,**f,

**res_S,**res_A,

norm_r=0,h, tau,

eps=0.1,

**A_u,**A_u_accur,

**A_u_u_accur,**r,

**z,**z_1d,*f_1d;

char fname[]="out.txt"; //имя выходного файла

setlocale (LC_CTYPE," Russian"); //ставим русский язык консоли

printf («Введите количество узлов сетки: «);//

scanf_s («%d» ,&n); //

n—; //

h=1./n; //вычисляем шаг сетки

printf («Вывод данных в файл %s…n», fname);

file = fopen (fname, «w»);

//инициализация переменных начало

mask=new int *[n+1];

u=new double *[n+1];

u_accur=new double *[n+1];

f=new double *[n+1];

res_S=new double *[n+1];

res_A=new double *[n+1];

A_u=new double *[n+1];

A_u_accur=new double *[n+1];

A_u_u_accur=new double *[n+1];

r=new double *[n+1];

z=new double *[n+1];

for (int i=0; i<=n; i++)

{

mask[i]=new int [n+1];

u[i]=new double [n+1];

u_accur[i]=new double [n+1];

f[i]=new double [n+1];

res_S[i]=new double [n+1];

res_A[i]=new double [n+1];

A_u[i]=new double [n+1];

A_u_accur[i]=new double [n+1];

A_u_u_accur[i]=new double [n+1];

r[i]=new double [n+1];

z[i]=new double [n+1];

}

for (int i=0; i<=n; i++)

{

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

{

mask[i][j]=0;

u[i][j]=0;

u_accur[i][j]=0;

f[i][j]=0;

res_S[i][j]=0;

res_A[i][j]=0;

z[i][j]=0;

}

}

//инициализация переменных конец

//создание маски квадратной области с вырезом и условием Дирихле на границах начало

int iStep, jStep;

iStep=n/2;

jStep=n/2;

for (int i=0; i<=n; i++)

{

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

{

if ((i>iStep) || ((i<=iStep)&&(j>jStep)))

{

mask[i][j]=ActualPoint;

}

if (i==0) mask[i][j]=FictivePoint; //левая граница

if (i==n) mask[i][j]=FictivePoint; //правая граница

if (j==0) mask[i][j]=FictivePoint; //нижняя граница

if (j==n) mask[i][j]=FictivePoint; //верхняя граница

}

}

//создание маски квадратной области с вырезом и условием Дирихле на границах конец

//вывод маски начало

fprintf (file," МАСКА ОБЛАСТИ! n");

for (int j=n; j>=0; j—)

{

for (int i=0; i<=n; i++)

{

fprintf (file," %d «, mask[i][j]);

}

fprintf (file," n");

}

//вывод маски конец

//вычисление u_accur, f в зависимости от маски начало

for (int i=0; i<=n; i++)

{

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

{

if (mask[i][j]==ActualPoint)

{

u_accur[i][j]=sin (M_PI*i*h)*sin (M_PI*j*h);

f[i][j] = -(-2*M_PI*M_PI*sin (M_PI*i*h)*sin (M_PI*j*h));

}

if (mask[i][j]==TopPoint)

{

u_accur[i][j] = sin (M_PI*i*h)*sin (M_PI*j*h);

f[i][j]=-(M_PI*sin (M_PI*i*h)*cos (M_PI*j*h)/h);

}

if (mask[i][j]==RightPoint)

{

u_accur[i][j] = sin (M_PI*i*h)*sin (M_PI*j*h);

f[i][j]=-(M_PI*cos (M_PI*i*h)*sin (M_PI*j*h)/h);

}

if (mask[i][j]==BottomPoint)

{

u_accur[i][j] = sin (M_PI*i*h)*sin (M_PI*j*h);

f[i][j]=(M_PI*sin (M_PI*i*h)*cos (M_PI*j*h)/h);

}

if (mask[i][j]==LeftPoint)

{

u_accur[i][j] = sin (M_PI*i*h)*sin (M_PI*j*h);

f[i][j]=(M_PI*cos (M_PI*i*h)*sin (M_PI*j*h)/h);

}

if (mask[i][j]==LeftBottomPoint)

{

u_accur[i][j] = sin (M_PI*i*h)*sin (M_PI*j*h);

f[i][j]=((M_PI*cos (M_PI*i*h)*sin (M_PI*j*h)/h)+(M_PI*cos (M_PI*j*h)*sin (M_PI*i*h)/h));

}

}

}

//вычисление u_accur, f в зависимости от маски конец

int matrixSize=0;

for (int i=0; i<=n; i++)

{

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

{

if (mask[i][j]≠FictivePoint)

{

matrixSize++;

}

}

}

z_1d=new double *[matrixSize];

for (int i=0; i<=matrixSize; i++)

{

z_1d[i]=new double[matrixSize];

}

//Проверка сопряженности — начало

int p=0,k=0;

for (int i=0; i<=n; i++)

{

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

{

if (i>0)

{

z[i-1][n]=0;

z[i-1][n-1]=0;

z[2][1]=0;

z[2][2]=0;

}

if (mask[i][j]≠FictivePoint)

{

z[i][j-1]=0;

z[i][j]=1;

A (n, mask, z, res_A);

for (int l=0; l<=n; l++)

{

for (int m=0; m<=n; m++)

{

if (mask[l][m]≠FictivePoint)

{

z_1d[p][k]=res_A[l][m];

p++;

}

}

}

k++;

p=0;

}

}

}

f_1d=new double [matrixSize];

p=0;

for (int l=0; l<=n; l++)

{

for (int m=0; m<=n; m++)

{

if (mask[l][m]≠FictivePoint)

{

f_1d[p]=f[l][m];

p++;

}

}

}

print (matrixSize, z_1d, f_1d);

//Проверка сопряженности — конец

//Проверка семетричночти и положительноопределенности оператора начало

Simetric (matrixSize, z_1d);

More_then₀(matrixSize, z_1d);

//Проверка семетричночти и положительноопределенности оператора конец

//приведение к треуг виду с правой частью начало

double d=0;

for (int k=0;k

{

for (int i=k+1;i

{

if (z_1d[k][k]==0)

{

SetConsoleTextAttribute (GetStdHandle (STD_OUTPUT_HANDLE), 4) ;

printf («Произошло деление на 0! nДля выхода нажмите любую клавишу… n»);

_getch ();

exit (0);

}

else

{

d = z_1d[i][k]/z_1d[k][k];

for (int j=0;j

{

z_1d[i][j]=z_1d[i][j]-z_1d[k][j]*d;

}

f_1d[i]=f_1d[i]-f_1d[k]*d;

}

}

}

print (matrixSize, z_1d, f_1d);

//приведение к треуг виду с правой частью конец

do

{

A (n, mask, u, res_A);

Substruction (n, mask, res_A, f, res_S);

A (n, mask, res_S, res_A);

tau=Scalar (n, mask, res_A, res_S)/Scalar (n, mask, res_A, res_A);

for (int i=0;i<=n;i++)

{

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

{

if (mask[i][j]≠0)

{

u[i][j]=u[i][j]-tau*res_S[i][j];

}

}

}

norm_r=Norm (n, mask, res_S);

count_it++;

}

while (norm_r > eps);

fprintf (file," Норма = %.40lf n", norm_r);

fprintf (file," nn");

fprintf (file," Количество итераций: %dn", count_it);

//проверка линейности оператора начало

A (n, mask, u, res_A);

copy (n, res_A, A_u);

//print (n, A_u);

A (n, mask, u_accur, res_A);

copy (n, res_A, A_u_accur);

//print (n, A_u_accur);

Substruction (n, mask, A_u, A_u_accur, res_S);

copy (n, res_S, r);

//print (n, res_S);

Substruction (n, mask, u, u_accur, res_S);

A (n, mask, res_S, res_A);

copy (n, res_A, A_u_u_accur);

//print (n, A_u_u_accur);

Substruction (n, mask, A_u_u_accur, r, res_S);

double maxDiff = 0;

for (int i=0;i<=n;i++)

{

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

{

if (abs (res_S[i][j])>maxDiff)

{

maxDiff = abs (res_S[i][j]);

}

}

}

fprintf (file," max (Lin) = %.25fn", maxDiff);

//проверка линейности оператора конец

fclose (file);

system («pause»);

return 0;

}

4. Результаты работы программы При n=6:

1) Внешний текстовый файл out. txt:

МАСКА ОБЛАСТИ!

0 0 0 0 0 0

0 1 1 1 1 0

0 1 1 1 1 0

0 0 0 1 1 0

0 0 0 1 1 0

0 0 0 0 0 0

100 -25 -25 -0 -0 -0 -0 -0 -0 -0 -0 -0 | 11

— 25 100 -0 -25 -0 -0 -0 -0 -0 -0 -0 -0 | 7

— 25 -0 100 -25 -0 -0 -25 -0 -0 -0 -0 -0 | 18

— 0 -25 -25 100 -0 -0 -0 -25 -0 -0 -0 -0 | 11

— 0 -0 -0 -0 100 -25 -0 -0 -25 -0 -0 -0 | 11

— 0 -0 -0 -0 -25 100 -25 -0 -0 -25 -0 -0 | 18

— 0 -0 -25 -0 -0 -25 100 -25 -0 -0 -25 -0 | 18

— 0 -0 -0 -25 -0 -0 -25 100 -0 -0 -0 -25 | 11

— 0 -0 -0 -0 -25 -0 -0 -0 100 -25 -0 -0 | 7

— 0 -0 -0 -0 -0 -25 -0 -0 -25 100 -25 -0 | 11

— 0 -0 -0 -0 -0 -0 -25 -0 -0 -25 100 -25 | 11

— 0 -0 -0 -0 -0 -0 -0 -25 -0 -0 -25 100 | 7

100 -25 -25 -0 -0 -0 -0 -0 -0 -0 -0 -0 | 11

0 94 -6 -25 0 -0 -0 -0 -0 -0 -0 -0 | 7

0 0 93 -27 0 -0 -25 -0 -0 -0 -0 -0 | 18

0 0 0 86 0 -0 -7 -25 -0 -0 -0 -0 | 11

0 0 0 0 100 -25 -0 -0 -25 -0 -0 -0 | 11

0 0 0 0 0 94 -25 -0 -6 -25 -0 -0 | 18

0 0 -0 0 0 0 86 -27 -2 -7 -25 -0 | 18

0 0 -0 0 0 0 -0 84 -1 -2 -8 -25 | 11

0 0 -0 0 0 0 -0 0 93 -27 -1 -0 | 7

0 0 -0 0 0 0 -0 0 0 85 -27 -1 | 11

0 0 -0 0 0 0 -0 0 0 -0 83 -28 | 11

0 0 -0 0 0 0 -0 0 0 -0 0 83 | 7

Норма = 0.81 934 754 114 870 502 430 768 034 217 984

Количество итераций: 65

max (Lin) = 0.71 054 273 576

2) Консоль

1. Демидович Б. П., Марон И. А., Шувалова Э. З. «Численные методы анализа. Приближение функций, дифференциальные и интегральные уравнения», М.: Наука, 1967. — 368 с.

2. Вержбицкий В. М. «Численные методы. Математический анализ и обыкновенные дифференциальные уравнения», М.: Высшая школа, 2001. — 384 с.

3. Бундаев В. В., Дамбаев Ж. Г. «Методические указания и контрольные задания по численным методам», Улан-Удэ: РИО ВСГТУ, 2003. — 16 с.