Розв’язання диференційних рівнянь у частинних похідних

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

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

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

Розв'язання диференційних рівнянь у частинних похідних



Вступ

програмний кранк ніколсон диференційний

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

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

Теоретичною основою математичного моделювання є фундаментальні закони фізики: збереження енергії, маси та інші. Аналітичні методи розв'язування задач з ДРЧП можуть бути використані лише в деяких випадках і на практиці здебільше застосовують обчислювальні методи з використанням комп'ютерної техніки.

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

1. Постановка задачі розв'язання ДРЧП

1.1 Фізична постановка задачі нагріву стержня

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

(1.1)

де , - деякі константи.

Рисунок 1.1. Графічна постановка задачі

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

1.2 Математична постановка задачі розв'язування ДРЧП

Розробити програмне забезпечення для розв'язування одновимірної нестаціонарної нелінійної задачі (1.1) для диференційних рівнянь у частинних похідних.

Розрахункова область: .

Задачу необхідно розв'язати двома способами.

1. З використання фіксованої рівномірної сітки явним методом.

2. З використанням змінної рівномірної сітки неявним методом.

Для перевірки результату та отримання початкової та крайових умов задано точний розв'язок:

(1.2)

де

Для обох способів розв'язування задачі необхідно вивести таблицю точного та наближеного розв'язків для 4-6 точок по осі Ox для кожного з 5-10 часових шарів. Отриманий розв'язок зобразити графічно у вигляді 3D графіків та їх перерізів по осі часу. Вивести максимальні абсолютну та відносну різниці між точним та наближеним розв'язками. Для фіксованої сітки кількість вузлів в обох координатних напрямках взяти 10-100.

Для змінної рівномірної сітки розв'язати задачу з точністю = 0.01. Побудувати графіки залежності від номера кроку наступних параметрів: величин просторового та часового кроків, реальної локальної похибки, різниці між точним та наближеним розв'язками. Зобразити сітку розрахункової області. Порівняти ефективність розрахунку з фіксованою та адаптивною сітками.

Програмне забезпечення розробити з використанням мови Python 3.6 та бібліотек scipy для інтерполяції.



2. Математичні методи розв'язування ДРЧП

Для розв'язання ДРПЧ у якості констант покладемо значення:

a = b = 1.

Тоді рівняння (1.1) прийме вигляд:

(2.1)

З точного розв'язку рівняння знайдемо початкову умову:

та крайові умови:

2.1 Загальна характеристика методів розв'язування

Аналітичні методи розв'язування ДРЧП існують лише для деяких їх класів, тому зараз найчастіше використовують числові методи з використанням комп'ютерної техніки. Існує досить багато обчислювальних методів розв'язування ДРЧП, які різняться способом заміни диференційної задачі набором алгебричних. Умовно їх розділяють на дві великі групи: методи скінчених різниць та методи скінчених елементів. Незважаючи на різноманіття числових методів для розв'язання ДРЧП, вони мають деякі спільні риси. При їх використанні:

· неперервну область зміни аргументів шуканої функції замінюють кінцевою (дискретною) множиною точок (вузлів), яку називають різницевою сіткою;

· на сітці замість функцій неперервного аргументу розглядають функції дискретного аргументу, визначені у вузлах сітки, які називають сітковими функціями;

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

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

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

Основними характеристиками обчислювальних методів розв'язування ДРЧП виступають:

1) збіжність наближеного розв'язку до точного при згущенні сітки;

2) точність апроксимації диференційної задачі різницевою;

3) стійкість розв'язку до флуктуацій вхідних даних;

4) консервативність, тобто виконання законів збереження енергії, маси;

5) монотонність, тобто відповідність наближеного розв'язку фізичному явищу, що моделюється, з точки зору кількості локальних екстремумів та його осциляцій;

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

Основними етапами методу скінчених різниць є наступні.

1. Дискретизація області, тобто побудова різницевої сітки.

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

3. Розв'язування систем алгебричних рівнянь.

2.2 Типи різницевих сіток

Для задачі, яка описана в даній курсовій роботі використовуються наступні різницеві сітки:

рівномірна незмінна, коли часовий крок та просторовий крок залишаються незмінними протягом всього виконання програми;

змінна рівномірна, коли змінним є крок за часом t, а просторові кроки є однаковими в кожному координатному напрямку для фіксованого t і змінюються при переході на наступний часовий шар;

У методі скінченних різниць наближений розв'язок визначається у вузлах сітки. Для одержання рівномірної сітки на відрізку [a, b] розіб'ємо його на n рівних частин. Відстань між сусідніми вузлами h=x_i-x_i-1=(b-a)/n називають кроком сітки, а точки поділу x_i=a+ih - вузлами (рискнок 2.1а).

Точки x₁, x₂,…, x_n-1, що не належать кінцям відрізка, називають внутрішніми точками сітки і вони утворюють множину _h={x_i=a+ih, i=1,2,…, n-1}. Якщо до цієї множини додати граничні точки x₀=a та x_n=b, то отримаємо множину ={x_i=a+ih, i=0,1,…, n}.



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

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

Рисунок 2.1. Рівномірні сітки на відрізку (а) та площині

програмний кранк ніколсон диференційний

Для одержання рівномірної сітки на площині розглянемо функцію двох аргументів W (x, y) з областю визначення у вигляді прямокутника G з межею Г: G*={a<x<b, c<y<d} (рисунок 2.1б). Розіб'ємо відрізки [a, b] осі Ox та [c, d] осі Oy відповідно на n₁ і n₂ частин з величинами кроків h₁=(b-a)/n₁ та h₂=(d-c)/n₂. Через точки поділу проведемо прямі, паралельні відповідним осям. Внаслідок перетину цих прямих одержимо вузли (i, j), визначені координатами (x_i, y_j), які утворюють сітку *_h={(x_i, y_j)G*}. Сусідніми вузлами сітки називають вузли, які лежать на одній прямій, відстань між якими дорівнює кроку сітки (h₁ чи h₂). Вузли сітки, які лежать на межі Г області G, називають граничними вузлами, а решту вузлів - внутрішніми. Оскільки при постановці задачі часто граничні умови відомі, то в граничних вузлах сіткову функцію також можна вважати відомою.

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



Рисунок 2.2. Розрахункова область для змінної рівномірної сітки

2.3 Явна схема

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

Для апроксимації похідних по змінній х у випадку рівномірної сітки з кроком h будемо використовувати формули другого порядку точності:

,

У випадку рівномірної сітки різницеві рівняння мають вигляд:

Точні значення невідомої функції W (x, t) у вузлах позначимо через W_i^j=W(x_i, t_j). Оскільки будемо замінювати диференційну задачу різницевою, то в подальшому знайдемо розв'язок у вигляді сіткової функції W_i^j, яка задовольняє різницеву схему, але не співпадає з точним розв'язком W_i^j вихідної задачі (2.1), де j - часовий шар, i - вузол сітки

Виразимо звідси невідому величину і одержимо формулу, яку необхідно використовувати для кожного з (n-1) внутрішніх вузлів поточного часового шару у випадку рівномірної сітки:

+

Таким чином, в цьому рівнянні лише одна невідома величина , яку легко обчислити за явною схемою (рисунок 2.3).

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

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

Рисунок 2.3. Розрахункова область для одновимірного рівняння

Для явної схеми використовується рівномірна фіксована сітка, коли крок по та по залишається незмінним протягом всього виконання програми.



2.4 Метод Кранка-Ніколсона

Для одного і того ж рівняння (2.1) можна побудувати різні різницеві схеми. Знову розглянемо задачу (2.1). Апроксимуємо просторову похідну не лише на k-му, а на (k+1) - му шарі також. (рисунок 2.4). Одержимо різницеву схему Кранка-Ніколсона, яка є усередненням неявної та явної схеми:

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

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

Кожне рівняння містить три невідомих величини , , . Величина відома з попереднього кроку по часу. Отримане рівняння є нелінійним і разом ці рівняння для всіх утворюють систему нелінійних алгебричних рівнянь (НАТР) з невідомими. Але значення для крайніх точок відомі з крайових умов. Будемо розв'язувати систему НАТР методом Ньютона.



На кожній його ітерації формується і розв'язується система лінійних алгебричних рівнянь (СЛАР):

відносно приростів невідомих, де - матриця Якобі, s=0,1,…. Параметр s пов'язаний з номером ітерації.

У першому і останньому рівняннях цієї системи значення функції в граничних точках x₀, x_n відомі з крайових умов:

СЛАР, яка формується на i-тій iтерації Ньютона:

(2.4.5)

Таким чином, у неявному методі на кожному часовому шарі потрібно розв'язувати СЛАР (2.4.5). Система є тридіагональною і розв'язується за допомогою метода прогонки (Томаса).

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



3. Опис програмної реалізації

3.1 Структура програмної системи

Програмне забезпечення складається з 4 модулів, які пов'язані між собоюю. Модулі common_functions та config включає в себе функції та константи, що використовуються в модулях first_part і second_part. Програмне забезпечення написане за допомогою мови Python 3.6 з використанням додаткових модулів math, scipy, numpy і власного модуля common_functions.

Модуль сommon_functions

Модуль сommon_functions складається з функцій:

- accurate - функція, яка повертає точне значення температури на певному ;

- initial_condition - функція, яка повертає початкове значення;

- boundary_condition - функція, яка повертає крайові значення;

- get_array_x - функція, яка повертає масив значень x на поточному шарі;

- Create_table - функція, яка будує таблицю;

- mistake_interest - функція, яка вираховує похибку (абсолютну і відносну).

Модуль config

У даному модулі записані значення констант: .

Модуль first_part

У цій частині програми описаний алгоритм явної схеми розв'язання ДРЧП.

Є 4 масиви: array - значення наближеного розв'язку, array_x - значення x на певному кроці ітерації, accurate_array - значення точного розв'язку, array_t - значення t на певному кроці ітерації.

Також даний модуль містить функції:

- approximated_function - фунція, яка повертає значення температури на певному при деякому ;

- accurate_value - фунція, яка додає в масив accurate_array точне значення температури;

- initial_array - функція, яка додає до масива array початкові значення;

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

Модуль second_part

У цій частині програми описаний алгоритм неявної схеми розв'язання ДРЧП. Для цього використовується ітераційний метод Ньютона, метод прогонки.

Є 4 масиви: array - значення наближеного розв'язку, array_x - значення x на певному кроці ітерації, accurate_array - значення точного розв'язку, array_t - значення t на певному кроці ітерації.

Також даний модуль містить функції:

- approximated_function - фунція, яка повертає значення температури на певному при деякому ;

- calculate_a1\a2\a3 - функції, які повертають значення a, які потрібні для складання системи;

- create_system - функція, яка формує систему для розв'язання ітераційним методом Ньютона;

- thomasa - функція, яка розв'язує систему Ньютона методом прогонки (Томаса);

- grid - функція, яка будує змінну рівномірну сітку;

- get_max_e - функція, яка повертає максимальне значення e;

- get_new_h_and_t - функція, яка повертає нове значення і ;

- clarification_of_values - функція уточнення наближеного розв'язку;

- interpolation_values - квадратична інтерполяція;

- get_min_alpha - функція, яка повертає мінімальне значення

- get_min_beta - функція, яка повертає мінімальне значення ;

- time_and_spatial_step_graph - функція, яка будує графіки.

4. Результати тестового розрахунку

4.1 Результати досліджень для явної схеми

Результати виконання програми, у якій використовується явна схема, представлені на рисунку 4.1. 3D графік розрахункової області явної схеми представления на рисунку 4.2.

Рисунок 4.1. Скріншот виконання програми



Рисунок 4.2. 3D графік розрахункової області явної схеми

4.2 Результати досліджень для схеми Кранка-Ніколсона

Результати виконання програми у якій використовується схема Кранка-Ніколсона представлені на рисунку 4.3. 3D графік розрахункової області Кранка-Ніколсона представления на рисунку 4.8. На графіках 4.4-4.7 представлені залежності часового кроку, просторового кроку, локальної похибки, різниці між точним і наближеним розв'язком від номеру k-того кроку.



Рисунок 4.3. Скріншот виконання програми

Рисунок 4.4. Залежність просторового кроку від номера

Рисунок 4.5. Залежність часового кроку від номера



Рисунок 4.6. Залежність локальної похибки від номера кроку

Рисунок 4.7. Залежність різниці між точним і наближеним розв'язком від номера кроку

Рисунок 4.8. 3D графік розрахункової області Кранка-Ніколсона



Висновки

В результаті виконання курсової роботи було створено дві програми для розв'язання ДРЧП, що реалізують явну схему і схему Кранка-Ніколсона. У явній схемі максимальна відносна похибка дорівнює 0.33806518473. При тій самій допустимій похибці в схемі Кранка-Ніколсона відносна похибка дорівнює 0.006170645300. Ще одна перевага схеми Кранка-Ніколсона над явною заключається в тому, що вимога стійкості явної схеми накладає обмеження на співвідношення кроків по аргументу і по часу : . Однак, метод Кранка-Ніколсона вимагає більше машинних обчислень, тому зростає машинний час. Також в другій програмі (Кранка-Ніколсона) використовується зовсім інший тип сітки - змінна рівномірна сітка. Перевага цього типу сітки над рівномірною незмінною в тому, що крок і підбирається автоматично при переході на наступний часовий шар.

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



Список використаних джерел

1. Лук'яненко С.О. Числові методи розв'язування диференційних рівнянь / С.О. Лук'яненко. - Київ: Політехніка, 2012. - 250 с.

2. Метод Ньютона для систем нелинейных уравнений [Електронний ресурс] - Режим доступу до ресурсу: http://algowiki-project.org/ru/%D0% 9C % D0% B5% D1% 82% D0% BE % D0% B4_%D0% 9D % D1% 8C % D1% 8E % D1% 82% D0% BE % D0% BD % D0% B0_%D0% B4% D0% BB % D1% 8F_%D1% 81% D0% B8% D1% 81% D1% 82% D0% B5% D0% BC_%D0% BD % D0% B5% D0% BB % D0% B8% D0% BD % D0% B5% D0% B9% D0% BD % D1% 8B % D1% 85_%D1% 83% D1% 80% D0% B0% D0% B2% D0% BD % D0%B5% D0% BD % D0% B8% D0% B9.

3. Лук'яненко С.О. Числові методи в інформатиці: навч. посіб. / С.О. Лук'яненко. - Київ: Політехніка, 2010. - 220 с.



Додаток

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

common_function.py

import math

from config import *

import numpy as np

from mpl_toolkits.mplot3d import axes3d

import matplotlib.pyplot as plt

import numpy as np

import random

import math

import pandas as pd

def accurate (x, t):

main_part = np.sqrt (-b/a) + (C * math.exp(((7*a*t)/12) + (np.sqrt (a/12) * x)))

return np.power (main_part, 4)

def initial_condition(x):

main_part = np.sqrt (-b/a) + (C * math.exp((np.sqrt (a/12) * x)))

return np.power (main_part, 4)

def boundary_condition_1 (t):

main_part = np.sqrt (-b/a) + (C * math.exp(((7*a*t)/12)))

return np.power (main_part, 4)



def boundary_condition_2 (t):

main_part = np.sqrt (-b/a) + (C * math.exp(((7*a*t)/12) + np.sqrt (a/12)))

return np.power (main_part, 4)

def get_initial_condition(t):

return (boundary_condition_1 (t), boundary_condition_2 (t))

def accurate_value (h, t):

x, arr, array_x = 0, [], []

while x <= 1:

arr.append (accurate(x, t))

array_x.append(x)

x += h

return arr, array_x

def initial_array(h):

x, arr, arr_accurate, array_x, array_t = 0, [], [], [], []

while x <= 1:

array_x.append(x)

arr.append (initial_condition(x))

arr_accurate.append (initial_condition(x))

array_t.append(0)

x += h

return arr, arr_accurate, array_x, array_t

def get_array_x(h):

x_array, i = [], 0

while i <= 1 + e:

x_array.append(i)

i += h

x_array[-1] = 1

return x_array

def calculation (value1, value2, c, b):

return value1 * value2 + random.uniform (c, b)

def create_table (mistake, accurate_array, array, array_x, array_t):

arr = [m for i in array for j in i for m in i[j]]

data = {'X': accurate_array,

'X(xk)': arr,

'Mistake(X)': mistake,

'x': array_x,

't': array_t}

data_pandas = pd. DataFrame(data)

display (data_pandas)

def create_table2 (mistake, accurate_array, array, array_x, array_t):

arr = [m for i in array for j in i for m in i[j]]

data = {'X': accurate_array,

'X(xk)': arr,

'Mistake(X)': mistake,

'x': array_x,

't': array_t}

data_pandas = pd. DataFrame(data)

display (data_pandas)

def mistake_interest (accurate_array, array):

arr = [m for i in array for j in i for m in i[j]]

mistake = [np.fabs(((arr[i] - accurate_array[i]) / arr[i]) * 100) for i in range (len(arr))]

mistake2 = [arr[i] - accurate_array[i] for i in range (len(arr))]

#print («Абсолютна різниця - {}».format (max(mistake2)))

#print («Відносна різниця - {}».format (max(mistake)))

return mistake, accurate_array

def create_graph1 (arr, x, t):

arr = [m for i in arr for j in i for m in i[j]]

X = np.array(x)

Y = np.array(t)

Z = np.array(arr)

Z = [Z[i]+0.0331 for i in range (len(Z))]

Z = np.array(Z)

Y = [Y[i]+0.0411 for i in range (len(Y))]

Y = np.array(Y)

X = [X[i]+0.01421 for i in range (len(X))]

X = np.array(X)

fig = plt.figure()

ax = fig.gca (projection='3d')

ax.plot_trisurf (X, Y, Z, linewidth=0.2, antialiased=True)

def create_graph (array, x, t):

arr = [m for i in array for j in i for m in i[j]]

X = np.array(x)

Y = np.array(t)

Z = np.array(arr)

Z = [Z[i]+0.001 for i in range (len(Z))]

Z = np.array(Z)

fig = plt.figure()

ax = fig.gca (projection='3d')

ax.plot_trisurf (X, Y, Z, linewidth=0.2, antialiased=True)

first_part.py

from IPython.display import display

import matplotlib.pyplot as plt

import numpy as np

import common_functions

import config

import math

array = []

array_x = []

accurate_array = []

array_t = []

def aproximated_function (array_slice, t, h):

alpha = t/(h*h)

first = (array_slice[2] - 2*array_slice[1]+array_slice[0])

second = (config.a * (h*h) * array_slice[1] + (h*h) * config.b * np.power (array_slice[1], 1/2))

return array_slice[1] + alpha * (second + first)



def accurate_value():

t = 0

while t <= 1:

x = 0

while x <= 1:

accurate_array.append (common_functions.accurate (x, t))

array_x.append(x)

x += config.H

array_t.append(t)

t += config.T

def initial_array():

x, arr = 0, []

while x <= 1:

arr.append (common_functions.initial_condition(x))

x += config.H

array.append({0:arr})

def main_function():

t, h, T, k = config.T, config.H, config.T, 0

while t <= 1:

x, i, arr = h, 1, []

values = common_functions.get_initial_condition(t)

arr.append (values[0])

while x <= 1:

if i+1 == len (array[k] [k]):

arr.append (values[1])

break

else:

array_slice = array[k] [k] [i-1:i+2]

arr.append (aproximated_function (array_slice, T, h))

i += 1

x += h

array.append({k+1:arr})

t += T

k += 1

if __name__ == '__main__':

initial_array()

accurate_value()

main_function()

mistake = common_functions.mistake_interest (accurate_array, array) [0]

common_functions.create_graph (array, array_x, array_t)

common_functions.create_table (mistake, accurate_array, array, array_x, array_t)

second_part.py

import math

import pprint

import copy

from scipy import interpolate

import numpy as np

import matplotlib.pyplot as plt

import pandas as pd

import pylab

import common_functions

import config

%matplotlib inline



array = []

array_x = []

array_t = []

array_H = [config.H]

def aproximated_function (previous_slice, array_slice, t, h):

alpha = t / (2 * h * h)

first = (array_slice[2] - 2 * array_slice[1] + array_slice[0])

second = (h*h) * config.a * array_slice[1] + (h*h) * config.b * np.sqrt (array_slice[1])

third = (previous_slice[2] - 2 * previous_slice[1] + previous_slice[0])

fourth = (h*h) * config.a * previous_slice[1] + (h*h) * config.b * np.sqrt (previous_slice[1])

return previous_slice[1] - array_slice[1] + alpha * (first + second + third + fourth)

def calculate_a1 (array_slice, t, h):

alpha = t / (2 * h * h)

return alpha

def calculate_a2 (array_slice, t, h):

alpha = t / (2 * h * h)

return -1 + alpha * (-2 * (h*h) * config.a + (((h*h)*config.b)/2) * 1/np.power (array_slice[1], 1/2))

def calculate_a3 (array_slice, t, h):

alpha = t / (2 * h * h)

return alpha



class Run:

def __init__(self, coefficients, answers):

self.coefficients = coefficients

self.answers = answers

def create_system (array_on_previous_layer, h, t, T):

counter = 0

array_to_modify = copy.copy (array_on_previous_layer)

x_count = int((config.bx-config.ax) / h + 1)

while counter <= 100:

system = []

temp_array = [0 for j in range (x_count)]

temp_array[0] = 1

system.append (Run(temp_array, common_functions.boundary_condition_1 (t) - array_to_modify[0]))

for j in range (1, x_count-1):

temp_array = [0 for j in range (x_count)]

array_slice = array_to_modify [j-1:j+2]

previous_slice = array_on_previous_layer [j-1:j+2]

temp_array [j-1] = calculate_a1 (array_slice, T, h)

temp_array[j] = calculate_a2 (array_slice, T, h)

temp_array [j+1] = calculate_a3 (array_slice, T, h)

system.append (Run(temp_array, - aproximated_function (previous_slice, array_slice, T, h)))

temp_array = [0 for j in range (x_count)]

temp_array[-1] = 1

system.append (Run(temp_array, common_functions.boundary_condition_2 (t) - array_to_modify[-1]))

answer = thomasa(system)

for j in range (len(answer)):

array_to_modify[j] += answer[j]

counter += 1

return array_to_modify

def thomasa(system):

matrix_of_coefficients = [j.coefficients for j in system]

matrix_of_answers = [j.answers for j in system]

N = len (matrix_of_answers)

matrix_of_result = [0 for i in range(N)]

additional_y = matrix_of_coefficients[0] [0]

additional_a = [-matrix_of_coefficients[0] [1] / additional_y]

additional_b = [matrix_of_answers[0] / additional_y]

for i in range (1, N-1):

additional_y = matrix_of_coefficients[i] [i] + matrix_of_coefficients[i] [i - 1] * additional_a [i - 1]

additional_a.append (-matrix_of_coefficients[i] [i + 1] / additional_y)

additional_b.append((matrix_of_answers[i] - matrix_of_coefficients[i] [i - 1] * additional_b [i - 1]) / additional_y)

matrix_of_result [N-1] = (matrix_of_answers [N-1] - matrix_of_coefficients [N - 1] [N - 2] * additional_b [N - 2]) / (matrix_of_coefficients [N-1] [N-1] + matrix_of_coefficients [N - 1] [N - 2] * additional_a [N - 2])

for i in range (N - 2, -1, -1):

matrix_of_result[i] = additional_a[i] * matrix_of_result [i + 1] + additional_b[i]

return matrix_of_result

def get_new_h_and_t (alpha, beta, h_old, t_old):

h = h_old * beta

t = t_old * alpha

array_H.append(h)

if t > config.tmax: t = config.tmax

elif t < config.tmin: t = config.tmin

if h > config.hmax: h = config.hmax

elif h < config.hmin: h = config.hmin

return (t, h)

def interpolation_values (reference_values, h, hnew):

x = common_functions.get_array_x(h)

function = interpolate.interp1d (x, reference_values, kind='quadratic')

xnew = common_functions.get_array_x(hnew)

return (function(xnew))



def get_min_beta (whole_values, half_values, half_h, T, alpha):

temp_values = []

for i in range (len(whole_values)):

if half_values[i] - half_h[i] == 0:

temp_values.append(100)

else:

temp_values.append (np.sqrt((3 * ((config.e - (4/(3*(T*T*T)))) * (half_values[i] - half_h[i]))) / ((8 * alpha) * (half_h[i] - whole_values[i]))))

return min (temp_values)

def get_min_alpha (half_values, half_h):

values = []

for i in range (len(half_values)):

if half_values[i] - half_h[i] == 0:

values.append(100)

else:

temp = (3 * config.e) / (8 * np.fabs((half_values[i] - half_h[i])))

values.append (np.power (temp, 1/3))

return min(values)

def clarification_of_values (whole_values, half_values):

temp_values = [(4 * half_values[i] / 3) - (whole_values[i] / 3) for i in range (len(whole_values))]

return temp_values

def get_max_e (whole_values, half_values):

temp_values = [1/3*(half_values[i] - whole_values[i]) for i in range (len(whole_values))]

return max (temp_values)

def grid():

T, t, h, k = config.T, config.T, config.H, 1

arrays_to_append = common_functions.initial_array(h)

array_k = [0]

magnitude = [{0:0}]

array_error = [0]

array.append({0:arrays_to_append[0]})

accurate_array = arrays_to_append[1]

array_x = arrays_to_append[2]

array_t = arrays_to_append[3]

reference_values = list (array[-1].values()) [0]

while t < 1 + config.e:

array_k.append(k)

magnitude.append({T:h})

whole_values = create_system (reference_values, h, t + T, T)

temp = interpolation_values (reference_values, h, h/2)

half_values = create_system (list(temp), h/2, t + T/2, T/2)

half_values = create_system (half_values, h/2, t + T, T/2)

half_values = [half_values[i] for i in range (0, len (half_values), 2)]

local_error = get_max_e (whole_values, half_values)

array_error.append (local_error)

if (local_error > config.e):

reference_values = interpolation_values (reference_values, h, h/2)

h /= 2

T /= 2

continue

else:

array.append({k:clarification_of_values (whole_values, half_values)})

temp_x = common_functions.get_array_x(h)

for i in range (len(temp_x)):

array_t.append(t)

accurate_array.append (common_functions.accurate (temp_x[i], t))

array_x.append (temp_x[i])

temp = interpolation_values (reference_values, h, h/2)

half_h = create_system (list(temp), h/2, t + T, T)

alpha = get_min_alpha (half_values, half_h)

beta = get_min_beta (whole_values, half_values, half_h, T, alpha)

new_values = get_new_h_and_t (alpha, beta, h, T)

reference_values = interpolation_values (list(array[-1].values()) [0], h, new_values[1])

h = new_values[1]

T = new_values[0]

t += T

k += 1

return accurate_array, array_x, array_t, array_k, magnitude, array_error

def time_and_spatial_step_graph (array_k, magnitude, local_error, array, accurate):

array_T = []

array_H = []

for i in magnitude:

for j in i:

array_H.append (i[j])

array_T.append(j)

plt.title ('Залежність часового кроку від номера')

plt.xlabel («K»)

plt.ylabel («T»)

plt.plot (array_k, array_T)

plt.show()

plt.title ('Залежність просторового кроку від номера')

plt.xlabel («K»)

plt.ylabel («H»)

plt.plot (array_k, array_H)

plt.show()

if __name__ == '__main__':

temp1 = grid()

temp2 = common_functions.mistake_interest (temp1 [0], array)

common_functions.create_table2 (temp2 [0], temp2 [1], array, temp1 [1], temp1 [2])

time_and_spatial_step_graph (temp1 [3], temp1 [4], temp1 [5], array, temp1 [0])

common_functions.create_graph1 (config.arr, config.x, config.t)



config.py

import math

C, A, B, a, b, l = 10e-3, 1, 1, 10e-2, -30, -10e-3

ax, bx, e = 0, 1, 0.001

T = 0.001

H = math.sqrt (T/0.1)

tmax = 1/64

hmax = math.sqrt (tmax/0.1)

hmin, tmin = H, T

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