Размещено на http://www.allbest.ru/

ЗМІСТ

• Вступ
• 1. Метод сіток
• 1.1 Постановка задачі
• 1.2 Розв'язання задачі Діріхле для рівняння Лапласа методом сіток у прямокутної області
• 2. Текст програми

ВСТУП

Значна кількість задач фізики і техніки приводять до диференціальних рівняння в приватних похідних (рівняння математичної фізики). Усталені процеси різної фізичної природи описуються рівняннями еліптичного типу. Точні розв'язки крайових задач для еліптичних рівнянь вдається отримати лише в окремих випадках. Тому ці завдання вирішують в основному наближено. Одним з найбільш універсальних та ефективних методів, отримали в даний час широке поширення для наближеного рішення рівнянь математичної фізики, є метод кінцевих різниць або метод сіток.

Метод сіток або, інакше, метод кінцевих різниць найбільш поширений для наближеного рішення диференціальних рівнянь в приватних похідних. Основна ідея методу сіток полягає в тому, що диференціальне рівняння, початкові і граничні умови замінюються системою звичайно-різницевих рівнянь алгебри, наближено представляють дану крайову задачу. Розглянемо застосування методу сіток до вирішення задач теплопровідності на прикладі двомірної задачі.

Для лінійних рівнянь в приватних похідних одним з очевидних прийомів є застосування методу сіток. Тоді рівняння в приватних похідних замінюється системою лінійних алгебраїчних рівнянь з числом невідомих, залежних від числа взятих точок. Основною перешкодою тут є недостатній обсяг пам'яті існуючих машин. Для нелінійних рівнянь в приватних похідних застосування методу сіток призводить, як правило, до непереборним труднощам.

Застосування методу сіток дозволяє звести диференціальну крайову задачу до системи нелінійних в загальному випадку алгебраїчних рівнянь щодо невідомих вузлових значень функцій.

Особливості в застосуванні методу сіток до рівнянь в приватних похідних викликані насамперед тим,що для рівнянь різних типів можуть бути коректно поставлені завдання різного виду з початковими чи граничними умовами. Метод сіток тісно пов'язаний з типом диференціального рівняння і для рівнянь різних типів потрібно обирати, взагалі кажучи, сітки різних видів.

Однак при застосуванні методу сіток до розглянутої електростатичної задачі виникають деякі принципові труднощі. З найбільшою повнотою було досліджено застосування методу сіток до рівнянь еліптичного типу.

Це не служить перешкодою до застосування методу сіток для областей з неідеальними кордонами. Однак для обґрунтування потрібні інші, більш тонкі теореми. Зокрема рівномірна метрика, якою ми користуємося у цій главі, стає мало природною.

1. МЕТОД СІТОК

Методи скінченних різниць, методи сіток -- чисельні методи розв'язку інтегро-диференціальних рівнянь алгебри, диференціального, інтегрального числення, основані на заміні диференціальних операторів різницевими операторами, інтегралів -- сумами, а функцій неперервного аргументу -- функціями дискретного аргументу. Така заміна приводить до системи, взагалі кажучи, нелінійних алгебраїчних рівнянь, які зрештою зводяться до лінійної системи деяким ітераційним методом. Якщо початкова задача має вигляд

,

де -- циліндрова область інтеграції -- межа області , -- її основа, -- шукана вектор-функція і -- задані вектор-функції, -- просторовий векторний аргумент, і -- оператори (не обов'язково обмежені), то найпростіша схема інтеграції початкового рівняння має вигляд:

,

Тут -- сіткова функція, що є розв'язком різницевого рівняння, -- різницеві оператори, залежні від параметрів , сітки, , -- сіткова область, що апроксимує деяким чином область , -- її границя і -- сіткові функції, що апроксимують функції і відповідно. Окремим випадком схеми (2) є схема з вагами, коли , -- ваговий коефіцієнт.

Схема (2) називається двошаровою, оскільки вона зв'язує між собою значення і різницевого розв'язку на двох тимчасових шарах ; можливі також і багатошарові схеми. Якщо оператор , де -- одиничний оператор, оборотний, то схема (2) може бути представлена у вигляді, що розв'язується де оператор називаєстья оператором кроку різницевої схеми і враховує крайові умови, а -- функція, залежна від і . Кажуть, що оператор , залежний від параметра , апроксимує (приблизно) оператор , якщо при . Тут -- деяке еталонне сімейство функцій, на якому перевіряється апроксимація (наприклад, сімейство достатньо гладких функцій).

Схема (2) називаєтсья коректною, або стійкою, якщо , де означає норму оператора в деякому просторі. (див. простір абстрактний у функціональному аналізі), яке може залежати від .

Схема (2) апроксимує рівняння (1), якщо і . Для лінійних систем рівнянь встановлені теореми збіжності, що стверджують, що збіжність різницевого рішення до розв'язку початкового рівняння виходить з апроксимації і коректності (стійкості) різницевої схеми.

Якщо властивості апроксимації, стійкості і збіжності мають місце лише при деякому співвідношенні між параметрами сітки , , де , то їх називають умовними. Якщо ж ці властивості справедливі при будь-якому співвідношенні між і , то їх називають абсолютними.

Схема (2) називається явною, якщо , і неявною, якщо .

Схеми, що абсолютно збіжні, існують лише в класі неявних схем. Як правило, при відповідному виборі параметрів схеми (наприклад вагових коефіцієнтів) неявні схеми є абсолютно стійкими, вони допускають скільки завгодно великий крок . Але перетворення оператора ускладнює алгоритм. У разі одновимірних задач неявні схеми реалізують методом факторизації; вони є достатньо економічними. Для багатовимірних задач неявні економічні схеми одержують за допомогою дробових кроків методу, який зводить багатовимірні задачі до послідовності одномірних або простіших задач.

Для розв'язання стаціонарних задач застосовують метод стаціонування (встановлення), в якому стаціонарний розв'язок розглядається як межа нестаціонарного розв'язку із стаціонарним (або що встановлюються) краєвими умовами. Відповідно до цього стаціонарну задачу вирішують ітераційним методом, аналогічним різницевому методу інтеграції (2). На відміну від нестаціонарного випадку, оператор для ітераційного процесу повинен бути сильно стійкий, тобто повинен задовольняти умові , . При розв'язку нелінійних задач, особливо в механіці суцільного середовища, застосовують комбінації схем інтеграції з ітераційними методами (т.з. ітерації по нелінійності).

Суть методу полягає в наступному. Область безперервної зміни аргументів, замінюється дискретним безліччю точок (вузлів), яке називається сіткою або гратами. Замість функції неперервного аргументу розглядаються функції дискретного аргументу, визначені у вузлах сітки і звані сітковими функціями. Похідні, що входять в диференціальне рівняння і граничні умови, замінюються різницевими похідними, при цьому крайова задача для диференціального рівняння замінюється системою лінійних або нелінійних алгебраїчних рівнянь (сіткових або різницевих рівнянь). Такі системи часто називають різницевими схемами. І ці схеми вирішуються відносно невідомою сіткової функції.

Застосування методу сіток до вирішення гіперболічного рівняння вимагає, щоб вибрані кроки сітки по осях х і у задовольняли певним співвідношенням.

Застосування методу сіток для дослідження стійкості просторово-періодичних рухів //Чисельні методи динаміки в'язкої рідини - Новосибірськ ІТПМ СВ АН СССP, 1979 - С.

Схема для моделювання рівняння. Розглянемо застосування методу сітки спрямованої дії для вирішення конкретних завдань розробки нафтових родовищ. Для визначення температурного поля при закачуванні гарячого агента в пласт потрібно інтегрування рівняння (2.1) або (2.3) при певних початкових і граничних умовах. Нехай початкова температура пласта дорівнює нулю, іншими словами, відлік температури ведеться від її початкового значення. На кордоні х 0 температура закачуваного агента приймається TK 100 З і передбачається, що на деякій відстані від нагнітальної галереї температура залишається рівною її початкового значення.

Аналіз питань, пов'язаних з застосуванням методу сіток для вирішення конкретних рівнянь, проведемо для хвильового рівняння, при вирішенні якого в достатній мірі проявляються властиві методу сіток особливості.

Електрична схема для моделювання рівняння. Слід зазначити, що при вирішенні практичних завдань із застосуванням методу сітки спрямованої дії необхідно з'єднати між собою точки з однаковим потенціалом. Отже, сітка спрямованої дії має більш універсальне призначення, ніж звичайна електрична сітка. Вона містить операційні підсилювачі аналогових обчислювальних машин для виконання різних математичних операцій.

У роботі В. Ф. Баклановской розглянуті теоретичні питання, пов'язані із застосуванням методу сіток для чисельного розв'язання рівнянь одновимірної нестаціонарної фільтрації газу.

Для лінійних рівнянь в приватних похідних одним з очевидних прийомів є застосування методу сіток. Тоді рівняння в приватних похідних замінюється системою лінійних алгебраїчних рівнянь з числом невідомих,залежних від числа взятих точок. Основною перешкодою тут є недостатній обсяг пам'яті існуючих машин. Для нелінійних рівнянь в приватних похідних застосування методу сіток призводить, як правило, до непереборним труднощам.

Згідно з цими рекомендаціями, розроблено метод розв'язання системи диференціальних рівнянь (3.1) і (3.4), суть якого полягає в застосуванні методу сіток для просторової координати.

Рішення відповідних пружно завдань про полях місцевих деформацій здійснено в основному або обчислювальними засобами, або експериментально із застосуванням методу сіток, чутливих покриттів або муару.

Розглянемо застосування методу сіток до вирішення внутрішньої задачі Діріхле для рівняння Лапласа - рівняння еліптичного типу.

Порівняємо числа арифметичних операцій, необхідних для вирішення розглянутих електростатичних задач як методом інтегральних рівнянь, так і методом сіток. При застосуванні методу сіток до вирішення польової задачі необхідно штучно обмежити безмежну область, зайняту полем (див. гл.

На підставі проведеного аналізу та накопиченого досвіду можна відзначити, що при розрахунках несталих неізотермічних режимів трубопровідного транспорту газів в одномірної постановці найбільш ефективний метод сіток, причому перевага дають неявні схеми. Крім того, застосування методу сіток доцільніше, якщо використовується схема інтегрування не нижче другого порядку.

Розрахунки газотранспортних систем для оперативного управління припускають використання одновимірних рівнянь енергії і теплопровідності. У цьому розділі зроблені висновки про доцільність застосування методу сіток з неявній схемою для спрощених задач, коли до співвідношень гідродинаміки і нестаціонарного, теплообміну для середовища додається одномірне рівняння теплопровідності для стінки трубопроводу. У цьому випадку доводиться вирішувати задачу відшукання функцій р, Т, Тст, М методом послідовних наближень.

Зміна величини етах характеризує поступове зниження із зменшенням температури граничної пластичної деформації. Розв'язання відповідних пружно-пластичних задач про полях місцевих деформацій здійснюється в основному обчислювальними засобами або експериментально із застосуванням методу сіток, тензо-чутливих покриттів або муару. З наведених даних випливає, що з пониженням температури спостерігається більш раннє падіння величини бк, ніж руйнівних напружень. Це пояснюється чутливістю розкриття тріщини до пониження температури, пов'язаної з граничною пластичною деформацією і полегшенням проростання тріщин.

1.1 Постановка задачі

Постановка задачі: необхідно знайти неперервну функціюu=u(x,y),

що задовольняє всередині деякої кінцевої областіD еліптичному рівнянню

Лапласа:

діріхле задача лаплас

а на границі Г цій області - умові:

u(x,y)/Г= f(x,y)

деf(x,y) - задана неперервна на Г функція.

Вважаємо, що Г- простий кусково-гладкий контур.

Задачу пошуку розв'язання еліптичного рівняння 2-го порядку, за умовою, що розв'язок на границі приймає наперед задане значення, називається задачею Діріхле або першою крайовою задачею.

Для розв'язання:

а) оберемо квадратну сітку з кроком h;

б) для укладання різницевого рівняння скористуємося формулою;

в) визначаємо шаблон сітки. Множина вузлів , значення сіткової функції uik яких входять в різницеве рівняння, називається шаблоном. В різницевий аналог рівняння :

наявні елементи, що розташовані справа, зліва, зверху та знизу від елементу uik (рис. 1), тобто, шаблон має форму хреста і носить назву шаблон-хрест;

Рисунок 1 - Шаблон „хрест” для еліптичного рівняння

г) перетворимо рівняння (2) до виду (3):

Отримали систему лінійних алгебраїчних рівнянь (3). Доведено, що система (3) завжди сумісна і має єдине рішення. Розв'язок системи можна отримати методом Гауса. При наявності великої кількості вузлів доцільно використовувати метод прогонки або ітераційні методи.

Специфічним ітераційним методом для розв'язку системи (3) є метод осереднення Лібмана, який є варіантом методу Гауса-Зейделя:

а) обираємо початкові наближення ;

б) послідовні наближення для внутрішніх вузлів сіточної області розраховуємо по формулі:

в) ітераційний процес продовжують, поки не стане виконуватися нерівність:

де е - мале наперед задане число, що характеризує точність розрахунків.

Доведено, що для будь-якого кроку h процес Лібмана збігається до точного розв'язку незалежно від вибору початкового наближення.

Погрішність наближеного рішення має порядок О(h2).

Для вибору початкового наближення часто користуються таким правилом: оскільки значення шуканої функції uik завжди обмежено зверху і знизу:

розумно обирати за середнє арифметичне величин u(x,y) у 4-х граничних точках, що розташовані на однієї горизонталі та вертикалі з точкою ui,k.

Доведено, що погрішність від заміни диференціального рівняння (1) різницевими рівняннями (3) має порядок O(h2) і оцінюється за допомогою нерівності:

Різницева схема (3) має властивості стійності та збіжності.

1.2 Розв'язання задачі Діріхле для рівняння Лапласа методом сіток у прямокутної області

У цьому підрозділі будуть розглянуті деталі розв'язання задачі Діріхле для рівняння Лапласа методом сіток.

Приклад

Нехай функція u(х,у) задовольняє всередині прямокутної області

D = {(х,у)/0 ? х ? 4; 0 ? у ? 3) рівняння Лапласа з граничними умовами:

u|х=0=0; u|у=0=0;

u|х=4= u|y=3 = ;

Необхідно розв'язати задачу крайову задачу Діріхле, обравши квадратну сітку з кроком h = 1

Розв'язання

Схема граничної області та сітки наведена на рис. 2. Оскільки в Scilab не передбачено використання нульового індексу елементів матриць, починаємо нумерацію індексів з номеру 1. На рис. 2 наведено нумерацію вузлів у відповідності з цим.

Розв'язання задачі проведемо методом Лібмана. Алгоритм розв'язання буде таким:

а) задаємо вхідні дані:

1) координати прямокутної області (нижнього лівого та крайнього правого кутів) ;

2) кількість кроків по кожній зі змінних;

3) eps - точність розв'язання системи рівнянь;

Рисунок 2 - Гранична область і крайові умови задачі Діріхле

б) створюємо 4 підпрограми-функції для розрахунків крайових умов. В разі переходу до задачі з іншими крайовими умовами у них необхідно вносити відповідні зміни;

в) задаємо початкові наближення для внутрішніх точок. За них обираємо середні значення точок крайових умов, що лежать на осях х і у, що проходять через точку, тобто:

;

г) реалізуємо ітераційний метод Лібмана, до тих пір, поки для крайньої з внутрішніх точок буде виконуватися умова збіжності:

|

д) друкуємо результати

Для розв'язання задачі укладаємо скрипт . У скрипті?

наведені коментарі, що дозволяють зясувати зміст кожної дії. Лістинг 1

//********РОЗВЯЗАННЯ ЗАДАЧІ ДІРІХЛЕ******** //*******Підпрограми-функції для крайових умов function y=Nach_Niz(x) y= 0 endfunction

function y=Nach_Verx(x) y = x л 2 / 16 endfunction function z=Nach_Lev(y) z = 0 endfunction

function z=Nach_Prav(y) z = y л 2 / 9 endfunction

//******Початок розрахунків*******

// ***Уведення координат кутових точок платівки x0=0, y0=0; //Лівий нижній кут x1=4, y1=3; // Правий верхній кут nx=4; //Кількість кроків по осі х ny=3; //Кількість кроків по осі у dx=(x1-x0)/nx;// Крок по осі Х dy=(y1-y0)/ny;// Крок по осі у

u=zeros(ny+1,nx+1);//визначення й обнуління матриці u //**** Формування крайових умов на границі сітки***** for i = 2:nx+1

x = x0 + (i-1) * dx; u(1, i) = Nach_Niz(x); u(ny+1, i) = Nach_Verx(x);

end

for i = 2:ny+1

y = y0 + (i-1) * dy; u(i, 1) = Nach_Lev(y); u(i,nx+1) = Nach_Prav(y);

end

for i=2:ny; for j=2:nx;

u(i,j)=(u(i,1)+u(i,nx+1)+u(1,j)+u(ny+1,j))/4;

end

end

//Метод Лібмана

d=1

while d>1e-5 a = u(ny, nx) for i = 2:ny for k = 2:nx

u(i, k) = (u(i - 1, k) + u(i + 1, k) + u(i, k - 1) +

u(i, k + 1)) / 4; end

end

d = abs(a - u(ny, nx)); end

Після закінчення розрахунків розв'язок - вектори у , х, матрицю u можна вивести на екран Результати наведені у табл. 1.

Таблиця 1 - Результат розв'язання прикладу

	х=0	х=1	х=2	х=3	х=4
у=0	0	0	0	0	0
у=і	0	0,0405074	0,0898696	0,1331149	0,1111111
у=2	0	0,0722196	0,1858920	0,3314878	0,4444444
у=3	0	0,0625	0,25	0,5625	1

Для більш точного розв'язання можна розділити область на більш дрібну сітку. Для наочного представлення результатів доцільно користуватися побудовою трьохвимірного графіка за допомогою комп'ютерних засобів.

2. ТЕКСТ ПРОГРАМИ

# include < stdio.h >

# include < fstream.h >

# include < conio.h >

# include < iostream.h >

# include < math.h >

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 END

int 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 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 <13; i + +) {printf ("% g", arrayX [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 <33; i + +) {printf ("% g", arrayY [i]);}

printf ("n");}// End of Regulation

void 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 <

cout <

}

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 + +) {if (arrayP [i] [j]. BelongsToDh == 1) {

f <

"" <

f.close ();

}

Размещено на Allbest.ru