close

Вход

Забыли?

вход по аккаунту

Барабанова Александра Сергеевна. Исследование задачи опускания бугра грунтовых вод в однородных слоях под действием силы тяжести (посредством использования различных квадратур)

код для вставки
2
СОДЕРЖАНИЕ
ВВЕДЕНИЕ…………………………………………………………….……….….3
ГЛАВА 1. ПОСТАНОВКА ЗАДАЧИ………………………...…….………..….7
1.1.
Основная система интегрального и дифференциального уравнений…..7
1.2.
Численная схема для основной системы посредством метода
дискретных особенностей и формулы прямоугольников ……..……......9
1.3.
Численная схема для основной системы посредством метода
дискретных особенностей и формул Гаусса ………………………........10
ГЛАВА 2. ИССЛЕДОВАНИЕ ВЛИЯНИЯ ЧИСЛЕННЫХ МЕТОДОВ ДЛЯ
РЕШЕНИЯ ПОСТАВЛЕННОЙ ЗАДАЧИ ………………….……...………...16
2.1
Решение интегралов различными численными методами…………....16
2.2.
Решение интегрального уравнения различными численными методами
…………………………………………………………………...…...…….18
2.3.
Решение интегрального уравнения квадратурными формулами Гаусса
при различном числе узловых точек.………………….………….....…..27
2.4.
Применение различных численных методов для решения
интегрального уравнения в задаче оседания бугра грунтовых вод…...37
2.4.1. Применение метода прямоугольников…..……………………...…38
2.4.2. Применение квадратурной формулы Гаусса …….………..……...41
2.4.3. Исследование влияния численных методов …..………………......46
Заключение.….…...………………………….…………..……...………..............50
Список литературы…….…….………….………………………...……….…….53
Приложения………………………….………….…………………….….……....55
3
ВВЕДЕНИЕ
Обоснование выбора темы исследования и ее актуальность
В связи с возрастающим использованием математических задач, особой
востребованностью пользуются задачи математической физики, которые
позволяют построить математические модели для описания происходящих
процессов
и
получить
решения
этих
задач.
Одним
из
наиболее
востребованных являются задачи, описывающие природные явления и
позволяющие провести анализ происходящих процессов.
При поливе, дожде, разливе водоемов и таянии снега образуются бугры
грунтовых вод, которые растекаются с течением времени, В результате этого
процесса в областях, смежных с бугром грунтовых вод, происходит
изменение уровня грунтовых вод. Последнее может приводить к негативным
последствиям, связанным с затоплением или заболачиванием территорий, и
требует построения математических моделей процесса изменения уровня
грунтовых вод для оценки и предотвращения нежелательного влияния.[1]
П. Я. Кочина предложила свести моделирование процесса изменения
уровня грунтовых вод к решению задачи об эволюции резкой границы Lt,
отделяющей область сухого грунта от области, насыщенной водой. Было
показано, что давление р в сухом грунте, а значит и на подвижной границе Lt,
должно быть приблизительно постоянным. Кроме того, было установлено,
что при описании движения жидкости в грунте можно пренебрегать
инерционными членами и использовать линейный закон фильтрации (вектор
скорости имеет потенциал). Тогда, согласно закону сохранения массы,
скорость жидкости в направлении нормали к границе Lt должна быть
пропорциональна нормальной скорости vn смещения границы [2].
В последствии П.Я. Кочина, Л.A. Галин и Э. Г. Кучеренко, а затем и
другие исследователи методами теории функции комплексного переменного
предложили точные и приближенные (в виде рядов) решения конкретных
4
практических задач [2].
Степень разработанности проблемы
В данной работе анализируется задача о совместном движении двух
жидкостей
в
постановке
Лейбензона-Маскета
(модель
―поршневого‖
вытеснения), причем внешняя жидкость рассматривается как фиктивная;
после построения основной системы интегрального и дифференциального
уравнений выбирается предельный случай, в котором вязкость и плотность
внешней жидкости стремятся к нулю.
В
математической
постановке
данной
задаче
используются
интегральные уравнения, которые необходимо численно решать выбранным
численным методом. До недавнего времени данная задачи решалась с
использованием квадратурных формул прямоугольников при численном
интегрировании и численном решении интегральных уравнений. Однако для
повышения точности решения предлагается использовать для численного
интегрирования квадратурные формулы Гаусса. В работе используются
квадратурные формулы Гаусса, и полученные результаты сравниваются с
уже реализованным методом прямоугольников для данной задачи.
Так
же
внимание
уделяется
выбору
современного
языка
программирования, способного в процессе счета обеспечить высокую
скорость и точность решения при работе с большим объемом информации.
Для увеличения точности решения приходится использовать матрицы
больших размерностей, что на практике существенно замедляет скорость
расчетов, либо вообще не удается получить искомое решение.
Исходя из рассуждений, был выбран высокоуровневый
язык
программирования общего назначения – Python [3]. NumPy – расширение
языка Python, добавляющее поддержку больших многомерных массивов и
матриц, вместе с большой библиотекой высокоуровневых математических
функций для операции с этими массивами [4].
5
Предмет исследования
Содержание
программного
обеспечения
для
математического
моделирования процесса оседания бугра грунтовых вод в однородных слоях
под действием силы тяжести при наличии включений различной природы и
возможности его применения для исследования задачи.
Объект исследования
Исследование влияния выбора метода решения системы интегральных
уравнений на основе реальной задаче о математическом моделировании
процесса оседания бугра грунтовых вод под действием силы тяжести.
Цель работы
Исследование задачи и разработка программного обеспечения для
задачи математического моделирования процесса оседания бугра грунтовых
вод в однородных слоях под действием силы тяжести, в которой и
используются различные методы решения интегральных уравнений и
исследование влияния выбора метода на точность решения.
Основные задачи исследования:
1. Изучить литературу по рассматриваемому вопросу и сделать
заключение о необходимости и эффективности решения задачи;
2. Изучение и повторение литературы по теории численных методов
для решения интегральных уравнений.
3. Для разработанной модели получить численную схему для
реализации задачи исследования влияния различных методов решения
интегральных уравнения;
4. Разработать программный комплекс для численного счета на одном
из современных языков программирования;
6
5. Исследовать
проанализировать
полученное
графики
и
приближенное
таблицы,
решение
задачи
демонстрирующие
и
поведение
различных параметров приближенного решения задачи в зависимости от
выбранного метода;
6. Сделать выводы по результатам исследований.
Теоретическое и практическое значение работы
Работа имеет теоретическое и практическое значение, т.к. полученные
теоретические результаты могут быть использованы для корректировки
модели и в дальнейшем выбора метода дающего наиболее точное решение, а
разработанный комплекс программ может применяться для дальнейших
исследований и решений задач прикладного характера, в которых
необходимо решать интегральные уравнения.
Структура работы
Во
введении
обосновывается
актуальность
решаемой
задачи,
анализируются труды по данной тематике, ставится цель и задачи
исследования, а также основные направления исследования.
В первой главе ставится физическая задача, которая будет исследована,
дается ее математическая формулировка, спроектирована численная схема.
Так же приводятся сведения о квадратурной формулы Гаусса и построение
его численной схемы для поставленной нами задачи.
Во
различных
второй
главе
интегральных
проводится
уравнений
исследование
в
точности
зависимости
от
решения
выбранного
численного метода. В дальнейшем идет исследование влияния выбора метода
для
поставленной
нами
задачи
в
первой
главе,
математического
моделирования опускания бугра грунтовых вод под действием силы тяжести.
В заключении делается вывод о полученных результатах исследования
и приводится список литературы.
7
ГЛАВА 1. ПОСТАНОВКА ЗАДАЧИ
1.1. Основная система интегрального и дифференциального
уравнений.
В неоднородном изотропном слое проводимости
область, занятую одной жидкостью с вязкостью
окруженную другой жидкостью вязкостью
рассмотрим
и плотностью
,
и плотностью Жидкости
разделены гладкой кривой Lt. Область D совместной фильтрации жидкостей
может быть ограничена непроницаемой прямой Lj , разделяющей грунт или
непроницаемые породы, или прямой LP , разделяющей грунт и свободную
жидкость.[6] (рис. 1.1)
Рисунок 1 - Постановка задачи
8
Фильтрацию
жидкостей
описывает
квазипотенциал
,
удовлетворяющий эллиптическому уравнению
(
)
(
)
(1.1)
и граничным условия
(
)
(1.2)
(1.3)
(1.4)
Здесь П - потенциал массовых сил,
давление на LP , и v = Кgrad<
- функция, определяющая
скорость фильтрации.
Полагаем, что подвижная граница Lt в любой момент времени t> 0
задана параметрически
контура
Lt
(
- параметр. Тогда перемещение
),где
с течением времени
t
описывается дифференциальным
на
(1.5)
уравнением
а начальное положение известно:
(
)
(1.6)
Отметим, что система (1.1), (1.2) и (1.5) при µ 1 → 0 и
→0иР =1
совпадает с математической постановкой Кочиной для задачи об опускании
бугра грунтовых вод [2, с. 547-548](рис. 1.2).
Все выражения (1.1)—(1.6) записаны в безразмерных величинах. Связь
между характерными величинами имеет вид
Систему (1.1)—(1.5) решаем методом нанесения на границу
квазипотенциала двойного слоя [7, 8]. В результате с учетом того, что
и
→0
→ 0, получаем систему интегрального и дифференциального уравнений:
,
где
-
,
-
(1.7)
потенциал невозмущенного течения, V =
— скорость
9
невозмущенного течения,
,
-( )
( ) (
∫
двойного слоя;
параметр
( )
)
(
)
оператор
-
квазипотенциала
- ядро; Ф 1 - квазипотенциал стока
с полным расходом, равным -1; П = у - потенциал силы тяжести; оператор
,
-( )
(
.
)
(
,
)
циркуляцией равной -1;
-( )
(
)
/
( )
∫
–
скорость
(
)
вихря
с
полной
- функция тока этого вихря; е х и е у - единичные
орты. Граничные условия (1.3) на LI или (1.4) на LP учитываются подбором
функций Ф 1 и
.
1.2. Численная схема для основной системы посредством
метода дискретных особенностей и формулы прямоугольников
Построим дискретный аналог основной системы интегрального и
дифференциального уравнений (1.7). Для этого в каждый момент времени
t j ,j =
, границу
представим в виде п — 1 приблизительно равных
отрезков. В результате получим множество точек *(
)
+.
Считаем, что на каждом из отрезков плотность квазипотенциала двойного
слоя постоянна.
Следуя [7, 8], построим дискретный аналог основной системы
интегрального и дифференциального уравнений (1.7). Получим систему
линейных алгебраических уравнений размерности
вычисления элементарных смещений
и
подвижной границы в каждый
момент времени tj:
∑
выражений для
∑
10
̅̅̅̅̅̅̅̅̅̅
̅̅̅̅
(
)
Отметим, что система (1.8) построена в предположении, что в каждый
момент времени tj отрезки, образующие подвижную границу L t . , имеют
одинаковую длину, поэтому в процессе численных расчетов после каждого
элементарного смещения границы производилось ее переразбиение на
равные по длине части.
1.3. Численная схема для основной системы посредством
метода дискретных особенностей и формул Гаусса
Рассмотрим неоднородное уравнение Фредгольма 2 рода [9]:
( )
(
∫
) ( )
Задача состоит в том, чтобы имея ядро (
функцию
( )
(1.9)
) и функцию ( ), найти
( ). При этом существование решения и его множественность
зависит от числа , называемого характеристическим числом (обратное ему
называется собственным). Переменные
фиксированный отрезок ,
и
пробегают здесь некоторый
-.
При решении интегрального уравнения пользуемся методом конечных
сумм [10], идея которого заключается в замене определенного интеграла
конечной суммой с помощью одной из квадратурных формул
∫
где
( )
∑
- абсцисса точек отрезка ,
( ),
-
квадратурной формулы, не зависящей от
(
(1.10)
), - коэффициенты
( ). Заменяя приближенный
интеграл в уравнении Фредгольма (1.9) по формуле (1.10) и полагая
будем иметь
,
11
∑
где
( )
(
(
),
( ).
)
Получаем системы линейных алгебраических уравнений относительно
значений
. Решив системы методом квадратурной формулы Гаусса, мы
получаем таблицу приближенных значений
в точках
. Это позволит
записать приближенное решение уравнения (1.9) в виде
( )
( )
∑
Пусть отрезок интегрирования ,
(
) .
- непрерывной функции ( )разбит
на n равных частей точками
( – шаг разбиения,
)⁄ ). Обозначим через ( ) сплайн-функцию, аппроксимирующую
(
подынтегральную функцию ( ).[10]
Пусть
,
на
каждой
)расположено узлов
-(
подынтегральная
функция
Предположим, что функция
многочленом
части
(
разбиения
в
),
принимает
которых
значения ( )(
).
( ) на каждой i-й части аппроксимируется
,
( )степени
-(
). При этом на многочлен
( ) накладывается 2 ограничения:
а) значения многочлена и подынтегральной функции равны в узлах
интерполяции, т.е ( )
),(
( )(
);
б) определенный интеграл от функции
на отрезке ,
( )
-
выражается через значения подынтегральной функции ( ) в узлах в виде
их линейных комбинаций, т.е
∫
( )
(
)
(
)
(
)
Квадратурные формулы Гаусса для выбранной степени
(1.11)
сплайны
12
будут определены, если из условий а) и б) удастся найти m неизвестных
коэффициентов и координаты m узлов
(
).
Задача решается одновременно для всех n частей разбиения отрезка
-, если выразить
,
Положим ̃ ( )
Тогда (
)
,
)через переменную
-(
/ ̃( )
.
̃ ( )(
- так:
,
̃
( ( ))
) и соотношение (1.11)
),(
перепишем в виде
∫ ̃( )
̃ ( )
̃ ( )
̃ (
)
(1.12)
Выведем квадратурную формулу Гаусса с тремя узлами (m = 3). Для
. Функция ̃ ( )
этого необходимо определить шесть величин:
- многочлен степени p, общий вид которого
̃( )
.
(1.13)
Подставив соотношение (1.13) в (1.12) и учитывая, что ̃ ( )
̃ ( )(
),
(
получим
тождество
относительно
коэффициентов
):
∫ (
)
̃ (
̃(
)
̃ (
)
)
Шесть неизвестных будут определены однозначно из системы шести
уравнений. В общем случае степень p аппроксимирующего многочлена
13
всегда является нечетным числом и связана с числом узлов m соотношением
p=2m-1. В частности, для трех узлов имеем многочлен пятой степени (p=5).
Множители при
в левой части тождества вычисляем так:
(
∫
)
{
Приравнивая подобные выражения в левой и правой частях тождества
при одинаковых коэффициентах
. Получим следующие шесть
уравнений:
̃
{
̃
̃
̃
̃
̃
̃
̃
̃
̃
{̃
̃
̃
̃
̃
̃
̃
̃
(1.15)
Эту систему позволяет упростить следующее свойство ее решения:
неизвестные (
) системы 2m уравнений вида
{
(где
(
) являются нулями многочлена Лежандра
)
нули
многочлена
принадлежат
интервалу
( )
(-1,1)
и
расположены симметрично относительно середины интервала.
В рассматриваемом частном случае m=3 и
Находим нули многочлена из уравнения
уравнения
√
( )
(
).
. Подставив корни
в (1.15), получим систему трех линейных
14
уравнений относительно переменных̃ ̃ ̃ :
̃
̃
̃
̃√
̃ √
̃
{
̃
Очевидное решение этой системы есть ̃
̃
⁄ ̃
⁄ .
Теперь подставим значения в соотношение (1.11):
∫
( )
(
где
(
)
(
√
)
(
))
√
Итак, квадратурная формула Гаусса с тремя узлами записывается в
виде
∫
( )
∫
( )
∑
(
∑∫
(
)
( )
(
)
(
))
Если подынтегральная функция имеет непрерывную производную
шестого порядка, то для оценки погрешности формулы Гаусса с тремя
15
узлами можно использовать неравенство
( )
|
( )
( )|
(
)
|
( )
( )|
При вычислении интеграла до достижения заданной точности
методом двойного пересчета условие окончания вычислений имеет вид
| ( )
( ⁄ )|
,
где k=2m (m – число узлов в квадратурной формуле Гаусса). Полагают,
что ( )
( ⁄ ) с точностью .
16
ГЛАВА 2. ИССЛЕДОВАНИЕ ВЛИЯНИЯ ЧИСЛЕННЫХ
МЕТОДОВ ДЛЯ РЕШЕНИЯ ПОСТАВЛЕННОЙ ЗАДАЧИ
Расчеты выполнены с помощью языка программирования Python [11].
Он
является
интерпретируемым
объектно-ориентированным
высокоорганизованным языком программирования. Он прост, вместе с тем
очень гибок и выразителен. Его стандартная библиотека включает в себя
большой объем полезных функций.
2.1. Решение интегралов различными численными методами
Вначале
рассмотрим
реализацию
решения
интеграла
методом
квадратурной формулой Гаусса:
∫
Программная
реализация
квадратурной
(2.1)
формулы
Гаусса
для
интеграла(2.1) с заданием количества разбиения отрезка интегрирования
представлена в Листинге 1.
В строке 6 мы объявляем функцию, которая вычисляет интеграл для
заданного подотрезка. 15 – 17 – узлы и коэффициенты квадратурной
формулы Гаусса. 9 – 11 находим узлы для заданного отрезка с помощью
формулы. Строка 18 – задаем разбиение. 22 – вызов функции реализующая
решение интеграла на подотрезке
17
1. #!/usr/bin/python
2. # -*- coding: utf-8 -*3. import numpy as np
4. import matplotlib.pyplot as plt
5. def f(x):
6. return 1/(1+x**2)
7. def function(a,b):
8.
n = 3, Res = 0.0, h = (b - a)/2
9.
mas_A = [ ((1**(k+1)-(-1)**(k+1))/(k+1)) for k in range(n)]
10.
mas_t = [np.sqrt(3./5),0,-np.sqrt(3./5)]
11.
mas_T = [ [mas_t[i]**k for i in range(n)]for k in range(n)]
12. A = np.linalg.solve(mas_T,mas_A)
13. for i in range(n):
14.
s = (b + a) / 2 + (b - a) / 2 * mas_t[i]
15.
Res = A[i] * f(s) + Res
16. return Res * h
17.a = 0., b = 1.
18.n = eval(raw_input('Разбиение n= '))
19.h = (b - a) / n
20.Res1 = 0.0, Res = 0.0
21.for i in range(n):
22.Res1 = function(a + i * h, a + (i + 1) * h)
23.print (Res)
24.Res = Res + Res1
25.print (Res)
26.print (np.pi / 4)
Листинг 1 - Программная реализация квадратурной формулы Гаусса для
интеграла(2.1)
18
При запуске программ получаем:
1. Метод Гаусса при количестве узловых точек 3 и при n =1
приближенное решение равно 0.785267034991
2. Метод Гаусса при количестве узловых точек 3 и при n =2
приближенное решение равно 0.785397805596
3. Методом Гаусса при количестве узловых точек 3 и при n =6
приближенное решение равно 0.785398163238
При точном решении 0.7853981633974483, мы видим, что чем на большее
количество подотрезков мы разбиваем заданный интервал интегрирования,
тем получаем более точный ответ.
2.2. Решение интегрального уравнения различными
численными методами
Рассмотрим решение уравнения Фредгольма 2 рода.
( )
Программная
∫
реализация
( )
обобщенной
(2.2)
формулы
трапеций
для
интегрального уравнения Фредгольма 2 рода (2.2) представлена в Листингах
2, 3, 4 [17], [18].
Далее больше не будет использоваться метод трапеций, так как он дает
очень большую погрешность.
19
1 #!/usr/bin/python
2 # -*- coding: utf-8 -*3
4 import numpy as np
5 import matplotlib.pyplot as plt
6 def f(x,s):
7
return x*np.exp(x*s)
8 def g(x):
9
return np.exp(x) #f(x)
10 def prov(i,j,A,x,x1):
11
if(i==j):
12
return 1+(A*f(x,x1))
13
else:
14
return (A*f(x,x1))
Листинг 2 – Объявление функций для решения интегрального
уравнения (2.2)
17 a=eval(raw_input('Введите [a: '))
18 b=eval(raw_input('Введите ,b] '))
19 n=eval(raw_input('Разбиение n= '))
20 h=(b-a)/n
21 A=[h/2,h]
22 X=[ a+i*h for i in range(n+1)]
23 L=[ [ 0 for k in range(n+1)] for k in range(n+1)]
24 R=[ 0 for k in range(n+1)]
25 v=[ 0. for k in range(n+1)]
26 u=[ 0. for k in range(n+1)]
Листинг 3 – Объявление переменных для интегрального уравнения
(2.2)
20
27 for i in range(n+1):
28
for j in range(n+1):
29
if((j==0)or(j==n)):
30
31
32
33
L[i][j]=prov(i,j,A[0],X[i],X[j])
else:
L[i][j]=prov(i,j,A[1],X[i],X[j])
R[i]=(g(X[i]))
34 print L
35 print R
36 Y=np.linalg.solve(L,R)
37 print Y
38 print X
39 for i in range(n+1):
40
for j in range(n+1):
41
if((j==0)or(j==n)):
42
43
44
45
v[i]=(A[0]*f(X[i],X[j])*Y[j])+v[i]
else:
v[i]=(A[1]*f(X[i],X[j])*Y[j])+v[i]
v[i]=g(X[i])-v[i]
46 for i in range(n+1):
47
u[i]=1
48 plt.plot(X,v,"green")
49 plt.plot(X,u,"blue")
50 plt.legend(["Trapecia","Tochnoe"])
51 q=[0.2,0.8,0.9,1.1]
52 plt.axis(q)
53 plt.show()
Листинг 4 – Решение интегрального уравнения (2.2) с помощью
метода Трапеции.
21
Программная
реализация
метода
Симпсона
для
интегрального
уравнения Фредгольма 2 рода (2.2) представлена в Листингах 2, 5, 6.
16 a=eval(raw_input('Введите [a: '))
17 b=eval(raw_input('Введите ,b] '))
18 n=eval(raw_input('Разбиение n= '))
19 h=(b-a)/n
20 A=[h/3,4*h/3,2*h/3]
21 X=[ a+i*h for i in range(n+1)]
22 L=[ [ 0 for k in range(n+1)] for k in range(n+1)]
23 R=[ 0 for k in range(n+1)]
24 v=[ 0. for k in range(n+1)]
25 u=[ 0. for k in range(n+1)]
26 for i in range(n+1):
27
for j in range(n+1):
28
if((j==0)or(j==n)):
29
30
31
32
33
34
L[i][j]=prov(i,j,A[0],X[i],X[j])
elif((j%2)==0):
L[i][j]=prov(i,j,A[2],X[i],X[j])
else:
L[i][j]=prov(i,j,A[1],X[i],X[j])
R[i]=(g(X[i]))
35 Y=np.linalg.solve(L,R)
Листинг 5 – Программная реализация решение интегрального
уравнения (2.2) с помощью метода Симпсона.
22
41 for i in range(n+1):
42
for j in range(n+1):
43
if((j==0)or(j==n)):
44
45
46
47
48
49
v[i]=(A[0]*f(X[i],X[j])*Y[j])+v[i]
elif((j%2)==0):
v[i]=(A[2]*f(X[i],X[j])*Y[j])+v[i]
else:
v[i]=(A[1]*f(X[i],X[j])*Y[j])+v[i]
v[i]=g(X[i])-v[i]
50 for i in range(n+1):
51
u[i]=1
52 plt.plot(X,v,"green")
53 plt.plot(X,u,"blue")
54 plt.show()
Листинг 6 – Программная реализация построения графика решения
интегрального уравнения (2.2) и точного решения.
В строке 41-49 находим приближенного решения интегрального
уравнения с учетом коэффициентов Симпсона.
Предварительные расчеты:
Для получения коэффициентов Гаусса и узловых точек, необходимо
провести
некоторые
предварительные
расчеты.
Вначале
найдем
коэффициенты и точки для случая n=2.
В теории уже описывалось, как находить нужные нам компоненты.
∫
(
)
{
23
Подставим k=0,1 в интеграл выше и получим систему:
∑
∑
{
Полином Лежандра 2 степени есть
( )
Приравниваем
( )=0, находим
√
. Для определения
имеем
систему
{
√
Отсюда
Далее
√
.
полученные
значения
используем
в
нашей
программе.
Впоследствии, как уже говорилось выше, нам не придется проводить эти
вычисления, так как имеется таблица с уже полученными значениями.
Программная реализация метода квадратурных форм Гаусса для
интегрального уравнения Фредгольма 2 рода (2.2) представлена в Листингах
2, 7, 8.
24
12 a=eval(raw_input('Введите [a: '))
13 b=eval(raw_input('Введите ,b] '))
14 n=eval(raw_input('Разбиение n= '))
15 print "Введите ti"
16 mas_T=[ eval(raw_input('')) for i in range(n)]
17 print "Введите Ai"
18 mas_A=[ eval(raw_input('')) for i in range(n)]
19 X=[((b-a)/2+(b-a)/2*mas_T[i]) for i in range(n)]
20 A=[ (mas_A[i]/2) for i in range(n)]
21 print X,A
22 L=[ [ 0 for k in range(n)] for k in range(n)]
23 R=[ 0 for k in range(n)]
24 v=[ 0. for k in range(n)]
25 u=[ 0. for k in range(n)]
26 for i in range(n):
27
for j in range(n):
28
if (i==j):
29
30
31
32
L[i][j]=1+(A[j]*f(X[i],X[j]))
else:
L[i][j]=(A[j]*f(X[i],X[j]))
R[i]=(g(X[i]))
33 Y=np.linalg.solve(L,R)
Листинг 7 – Программная реализация решения интегрального уравнения (2.2)
с использованием квадратурной формулы Гаусса.
25
41 for i in range(n):
42
u[i]=1
43 plt.plot(X,v,"green")
44 plt.plot(X,u,"blue")
45 plt.show()
Листинг 8 – Программная реализация построения графика решения
интегрального уравнения (2.2) и точного решения.
В строке 33 находятся приближенные значения функции ( ).
При запуске программ получаем[12]:
Рисунок 2 – График решения интегрального уравнения (2.2) ,
полученного при использовании квадратурной формулы Гаусса (n=2)
26
Рисунок 3. График решения интегрального уравнения (2.2) ,
полученного при использовании обобщенной формулы Симпсона (n=2)
Рисунок 4 - График решения интегрального уравнения (2.2) ,
полученного при использовании обобщенной формулы трапеции (n=2)
27
Рисунок - 5. График решения интегрального уравнения (2.2) ,
полученного при использовании квадратурной формулы Гаусса и
обобщенной формулы Симпсона (n=6)
На рисунках 2, 3, 4 видим, как рассматриваемые нами функции ведут
себя относительно точного решения. На рисунке 5 идет сравнение 2 методов
относительно точного решения. На основе этих рисунков приходим к
выводу, что более точное решение дает квадратурная формула Гаусса.
2.3. Решение интегрального уравнения квадратурными
формулами Гаусса при различном числе узловых точек
Рассмотрим теперь решение уравнения Фредгольма 2 рода.
( )
∫
( )
(2.3)
Точное решение этого интегрального уравнения равно единице.
28
Программная
реализация
квадратурной
формулы
Гаусса
для
интегрального уравнения Фредгольма 2 рода (2.3) представлена в Листингах
9, 10, 11 [17], [18].
1 #!/usr/bin/python
2 # -*- coding: utf-8 -*3
4 import numpy as np
5 import matplotlib.pyplot as plt
6
7
8 font = {'family': 'Verdana',
9
'weight': 'normal'}
10 plt.rc('font', **font)
11
12 def f(x,s):
13 return x*np.exp(x*s)
14 def g(x):
15 return np.exp(x)
Листинг 9 – Программная реализация решения интегрального
уравнения(2.3), определения функций.
29
17 a = eval(raw_input('Введите [a: '))
18 b = eval(raw_input('Введите ,b] '))
19 n =eval(raw_input('Разбиение n = '))
20 s2 = 0., m = 2
21 c = [-0.57735,0.57735], A = [1,1]
22 L = [ [ 0 for k in range(m*n)] for k in range(m*n)]
23 R = [ 0 for k in range(m*n)]
24 Y = [ 0 for k in range(m*n)]
25 X = [ 0. for j in range(m*n)]
26 h = (b - a)/n; h1 = h/2; c1 = [c[i] * h1 for i in range(m)]
27 x2 = a + h1;
28 s2 = 0
29 k = 0
30 for i in range(n):
31 for j in range(m):
32
X[k] = x2 + c[j]*h1+(i*h)
33
k=k+1
34
35 for k in range(m*n):
36 for j in range(m*n):
37
38
39
40
if(k==j):
L[k][j] = 1. + A[0] * f(X[k],X[j]) * h1
else:
L[k][j] = A[0] * f(X[k],X[j]) * h1
41 R[k] = g(X[k])
42 Y = np.linalg.solve(L,R)
Листинг 10 – Программная реализация решения интегрального уравнения
(2.3), применения численного метода.
30
49 u = [ 0. for k in range(n*m)]
50 for i in range(n*m):
51 u[i] = 1.
52
53 plt.plot(X,Y,"green")
54 plt.plot(X,u,"blue")
55 q = [0.01,0.99,0.9999,1.0001]# для того, что бы можно было разглядеть
погрешность
56 plt.axis(q)
57 plt.title(u'При '+'n = '+str(n))
58 plt.legend([u"Гаусса",u"Точное"])
59 plt.show()
Листинг 11 – Программная реализация решение интегрального уравнения
(2.3), построения графика решения.
В итоге мы будем получаем m*n точек, поэтому в строках 26-28 и 49
задаем соответствующую погрешность. 30-33 вычисляются все х после
разбиения. 35–41 алгоритм в котором мы вычисляем левую и правую части и
находим неизвестные значения. 53-54 строим точное решение и
приближенное.
При запуске программ получаем[12]:
31
Рисунок 6 - График решения интегрального уравнения (2.3) ,
полученного при использовании квадратурной формулы Гаусса при
количестве узловых точек равном 2 и n = 1.
Рисунок 7 - График решения интегрального уравнения (2.3) ,
полученного при использовании квадратурной формулы Гаусса при
количестве узловых точек равном 2 и n = 2.
32
Рисунок 8 - График решения интегрального уравнения (2.3) ,
полученного при использовании квадратурной формулы при
количестве узловых точек равном 2 и n = 10. (точное и приближенное
решение сливаются)
На рисунках 6, 7, 8 видим, как рассматриваемая нами функции ведет
себя относительно точного решения. По этому примеру можно подвести
итог, что чем больше разбиение отрезка, тем точнее приближенное решение.
Рассмотрим следующий пример:
( )
Программная
∫
( )
реализация
(
квадратурной
)
формулы
(2.4)
Гаусса
для
интегрального уравнения Фредгольма 2 рода (2.4) представлена в Листинг
12, 13, 14[17], [18].
33
1 #!/usr/bin/python
2 # -*- coding: utf-8 -*3
4 import numpy as np
5 import matplotlib.pyplot as plt
6
7
8 font = {'family': 'Verdana',
9
'weight': 'normal'}
10 plt.rc('font', **font)
11
12 def f(x,s):
13 return (-1./2.)*np.exp(x*s)
14 def g(x):
15 return (1.-(1./(2.*x))*(np.exp(x)-1.))
Листинг 12 – Программная реализация решения интегрального
уравнения (2.4), определения функций.
34
17 a = eval(raw_input('Введите [a: '))
18 b = eval(raw_input('Введите ,b] '))
19 n = eval(raw_input('Разбиение n = '))
20 s2 = 0., m = 2
21 c = [-0.57735,0.57735], A = [1,1]
22
23 L = [ [ 0 for k in range(m*n)] for k in range(m*n)]
24 R = [ 0 for k in range(m*n)]
25 Y = [ 0 for k in range(m*n)]
26
27 X = [ 0. for j in range(m*n)]
28 h = (b - a)/n; h1 = h/2; c1 = [c[i] * h1 for i in range(m)]
29 x2 = a + h1;
30 s2 = 0, k = 0
31 for i in range(n):
32 for j in range(m):
33
X[k] = x2 + c[j]*h1+(i*h)
34
k=k+1
35 for k in range(m*n):
36 for j in range(m*n):
37
38
39
40
if(k==j):
L[k][j] = 1. + A[0] * f(X[k],X[j]) * h1
else:
L[k][j] = A[0] * f(X[k],X[j]) * h1
41 R[k] = g(X[k])
42 Y = np.linalg.solve(L,R)
Листинг 13 – Программная реализация решения интегрального
уравнения (2.4), применение квадратурной формулы Гаусса.
35
49 u = [ 0. for k in range(n*m)]
50 for i in range(n*m):
51 u[i] = 1.
52
53 plt.plot(X,Y,"green")
54 plt.plot(X,u,"blue")
55 q = [0.01,0.99,0.9999,1.0001]# для того, что бы можно было разглядеть
погрешность
56 plt.axis(q)
57 plt.title(u'При '+'n = '+str(n))
58 plt.legend([u"Гаусса",u"Точное"])
59 plt.show()
Листинг 14 – Программная реализация решения интегрального уравнения
(2.4), построение графика решения.
При запуске программ получаем:
Рисунок 9 - График решения интегрального уравнения (2.4) , полученного
при использовании квадратурной формулы Гаусса при количестве узловых
точек равном 2 и n = 1.
36
Рисунок 10 - График решения интегрального уравнения (2.4) ,
полученного при использовании квадратурной формулы Гаусса при
количестве узловых точек равном 2 и n = 2.
Рисунок 11 - График решения интегрального уравнения (2.4) ,
полученного при использовании квадратурной формулы при
количестве узловых точек равном 2 и n = 10. (точное и приближенное
решение сливаются)
37
На рисунках 9, 10, 11 видим, как рассматриваемая нами функции ведет
себя относительно точного решения. По этому примеру можно подвести
итог, что чем больше разбиение отрезка, тем точнее приближенное решение.
2.4. Применение различных численных методов для решения
интегрального уравнения в задаче оседания бугра грунтовых
вод
В [2] приводится исследование плоской задачи о растекании бугра
грунтовых вод методом решения нелинейного уравнения на границе. Так, для
бугра, в начальный момент времени представленного параметрическими
уравнениями [2, с. 554]
при t = 0
(2.5)
и расположенного в безграничном однородном грунте, закон опускания
частицы контура, находящейся на вершине бугра, имеет вид
(
)
.
(
√
)
(
) /.
(2.6)
Для численного решения задачи используется система (1.8). В этой
системе плоскому случаю однородного грунта соответствуют значение
проводимости
(
)
P,
равное
1,
и
функции
(
)
и
.
В статье [1] проведено сопоставление результатов численного решения
задачи об оседании бугра (2.6) с точным решением (2.5). Результаты
представлены в табл. 1. статьи [1].
Решение (2.6) позволяет протестировать процесс изменения положения
границы раздела воды и воздуха
.
Для тестирования результатов
вычисления поля скоростей по схеме (1.8) в статье воспользовались
приближенным решением Кочиной [8, с. 559]:
38
(
Решение
(2.7)
)
описывает
.
эволюцию
(2.7)
бугра
грунтовых
вод,
подчиняющуюся следующему закону:
.
/
.
Всеми этими выводами мы будем пользоватесь при реализации
программного обеспечения для исследования влияния различных методов
решения интегральных уравнений на данную задачу.
2.4.1. Применение метода прямоугольников
Ниже приведен основной код программы[11,12]
1.
thet=np.linspace(-20,20,n+1,endpoint=True)
2.
x=np.array([thet,H_K/np.pi*(np.arctan((thet+R_K)/(c_K*t_K))np.arctan((thet-R_K)/(c_K*t_K)))+d2])
3.
plot(x[0],x[1])
4.
isMove=True; T=0; cycle=0
Листинг 15 – Программная реализация решения задачи опускания
бугра грунтовых вод с использованием метода прямоугольников (задание
начальных узлов).
39
5.
while isMove:
6.
A,f=ieFredgolm2dk(x,x,lambdat,alphat,q,z,varphiSo
urces,Py,Phi1)
7.
g=np.linalg.solve(A,f)
8.
9.
W=VSources(q,x,z,Kx,Phi1)+VDL(x,x,g,Hx,Psi2)
10.
if cycle==0: quiver(x[0],x[1],W[0],W[1])
11.
x+=W*dt
12.
13.
x[0][0]=-20; x[0][n]=20
14.
interpolant=interp1d(x[0].ravel(), x[1].ravel())
15.
x1new=interpolant(thet)
16.
x=np.array([thet,x1new])
17.
T+=dt; cycle+=1;
18.
if cycle%50==0: print (cycle) ; plot(x[0],x[1])
19.
ymax=np.max(x[1])
20.
print ("%12.10f" % T, ymax)
21.
eps_H=0.00000000001
22.
ymaxH=np.abs(ymax-H_min)
23.
if (ymaxH<=eps_H or ymax<=H_min ): isMove=False; plot(x[0],x[1])
Листинг 16 – Программная реализация решения задачи опускания
бугра грунтовых вод с использованием метода прямоугольника.
Строки 1-2 задают начальные точки. В 5 строке начинается цикл в
котором происходит моделирование процесса опускания бугра. Решение
интегрального уравнения происходит в 6 строке, здесь мы переходим в
модуль, в котором написана реализация метода прямоугольника. 7 строка –
нахождение узлов. 9 строка – решение дифференциального уравнения, с
помощью которого задается смещение границы. Интерполяция – 14-15
строки и задание нового контура после смещения, 16 строка. Цикл
40
останавливается после того, как наша граница станет близка кграницы, ниже
которой она не может опуститься.
1. def
ieFredgolm2dk(xrw,xcl,lmd,alph,q,z,vphFun,Py=_Py_,Phi1=_Phi1_,isDiag
=1,minDist=myeps):
2.
nrw=xrw[0].shape[0]-1
3.
ncl=xcl[0].shape[0]-1
4.
x=np.array([(xrw[0,:nrw]+xrw[0,1:])/2,(xrw[1,:nrw]+xrw[1,1:])/2])
5.
y=np.array([(xcl[0,:ncl]+xcl[0,1:])/2,(xcl[1,:ncl]+xcl[1,1:])/2])
6. YC0,XC0=np.meshgrid(y[0,:],x[0,:])
7.
YC1,XC1=np.meshgrid(y[1,:],x[1,:])
8. Nrml=Normal(xcl)
9.
KF0=withoutZero(Py*sp.diff(Phi1,y0),minDist)(XC0,XC1,YC0,YC1)
10.
KF1=withoutZero(Py*sp.diff(Phi1,y1),minDist)(XC0,XC1,YC0,YC1)
11.
MtrA=KF0*Nrml[0]+KF1*Nrml[1]
12.
MtrA=-2*lmd*MtrA + ( np.eye(nrw) if isDiag else 0 )
13.
vtrf=2.*lmd*vphFun(q,x,z,Phi1)-2*alph*x[1]
return MtrA,vtrf
Листинг 17 – Модуль реализующий решение интегрального (метод
прямоугольника) и дифференциального уравнения.
Строка 4 – 5 определяет узлы по формуле прямоугольников. 9 – 10
подстановка этих узлов в наше уравнение. 12 – находим левую часть
уравнения, 13 – правую[12].
41
Рисунок 12 - График решения задачи математического моделирования
об опускании бугра грунтовых вод под действием силы тяжести при
800 точках, где интегральное уравнение решается с помощью формулы
прямоугольника.
На (рис. 12) мы видим, как с течением времени смещается кривая L. Но
по сложно сделать какие-то конкретные выводы, поэтому необходимо
заполнить таблицу получаемыми данными в ходе варьирования параметров
для этой задачи.
2.4.2. Применение квадратурной формулы Гаусса
Ниже приведен основной код программы[13,14]
42
1.
n=100, dt=0.1/4.0
2.
t_K=0.2; H_K=1.0;
3.
R_K=1.0; c_K=1.0
4.
m=8
5.
mas_T = [-0.96028986, -0.79666648, -0.52553242, -0.18343464,
0.18343464, 0.52553242, 0.79666648, 0.96028986]
6.
mas_A = [0.1012854, 0.22238104, 0.31370664, 0.36268378,
0.36268378, 0.31370664, 0.22238104,0.1012854]
7.
a = -20.; b = 20.
8.
h = (b - a)/n; h1 = h/2.
9.
x2 = a + h1;
10.
k = 0; end=a
11.
thet = np.arange(n*m)*1.
12.
Akoef = np.arange(n*m)*1.
13.
for i in range(n):
14.
end=end+h
15.
for j in range(m):
16.
thet[k] = x2 + mas_T[j]*h1 + (i*h)
17.
Akoef[k]=(end-(end-h))*(mas_A[j]/2)
18.
k=k+1
19.
20. x=np.array([thet,H_K/np.pi*(np.arctan((thet+R_K)/(c_K*t_K))np.arctan((thet-R_K)/(c_K*t_K)))+d2])
21. H0=H_K/np.pi*(np.arctan((0.+R_K)/(c_K*t_K))-np.arctan((0.R_K)/(c_K*t_K)))+d2
22.
r_skv = H0/5.0
23.
H_min = H0/5.0+d2
24.
plot(x[0],x[1])
Листинг 18 – Программная реализация решения задачи опускания
бугра грунтовых вод с использованием квадратурной формулы Гаусса,
объявление переменных, задание начальных узлов.
43
31.while isMove:
32.
A,f=ieFredgolm2dk(x,x,Akoef,h1,lambdat,alphat,q,z,varphiSources,Py,Phi1)
33.
g=np.linalg.solve(A,f)
34.
35.
W=VSources(q,x,z,Kx,Phi1)+VDL(x,x,g,Hx,Psi2)
36.
# if cycle==0: quiver(x[0],x[1],W[0],W[1])
37.
38.
x_new=x + W*dt
39.
x_new[0][0]=x[0][0]; x_new[0][n*m-1]=x[0][n*m-1]
40.
interpolant=interp1d(x_new[0].ravel(), x_new[1].ravel())
41.
x1new=interpolant(x[0])
42.
x=np.array([x[0],x1new])
43.
44.
T+=dt; cycle+=1;
45.
if cycle%50==0: print (cycle) ; plot(x[0],x[1])
46.
47.
ymax=np.max(x[1])
48.
print ("%12.10f" % T, ymax)
49.
eps_H=0.00000000001
50.
ymaxH=np.abs(ymax-H_min)
51.
52.
if cycle==500: isMove=False; plot(x[0],x[1])
53.
if (ymaxH<=eps_H or ymax<=H_min ): isMove=False; plot(x[0],x[1])
Листинг 19 – Программная реализация решения задачи опускания
бугра грунтовых вод с использованием квадратурной формулы Гаусса.
Строки 7-8 задают узлы и коэффициенты для квадратурной формулы
Гаусса. 9 – задание интервала с которым мы работаем. 10 – разбиваем наш
отрезок на n частей и для каждого отрезка находим значение который нужен
44
для получения узлов и коэффициентов для определенного нами интервала. 15
– 20 строка задают начальные точки применяя преобразования. Строка 22 –
подставляем полученные нами значения в формулу, которая задает наш
график. 23 – 25 задается Н_min с которым далее будет идти сравнение в
цикле. Задает минимальное значение y при котором имеет место решение. В
31 строке начинается цикл в котором происходит моделирование процесса
опускания бугра. Решение интегрального уравнения происходит в 32 строке,
здесь мы переходим в модуль, в котором написана реализация метода
прямоугольника. 33 строка – нахождение узлов. 35 строка – решение
дифференциального уравнения, с помощью которого задается смещение
границы. Интерполяция – 40-41 строки и задание нового контура после
смещения, 42 строка. Цикл останавливается после того, как наша граница
станет близка к границе, ниже которой она не может опуститься.
Ниже приведен модуль, отвечающий за решение интегрального
уравнения.
45
1.
def
ieFredgolm2dk(xrw,xcl,Akoef,h,lmd,alph,q,z,vphFun,Py=_Py_,Phi1=_Phi1_,isDi
ag=1,minDist=myeps):
2.
nrw=xrw[0].shape[0]-1
3.
ncl=xcl[0].shape[0]-1
4.
x=np.array([(xrw[0,:nrw]),(xrw[1,:nrw])])
5.
y=np.array([(xcl[0,:ncl]),(xcl[1,:ncl])])
6.
YC0,XC0=np.meshgrid(y[0,:],x[0,:])
7.
YC1,XC1=np.meshgrid(y[1,:],x[1,:])
8.
Nrml=Normal(xcl)
9.
KF0=withoutZero(Py*sp.diff(Phi1,y0),minDist)(XC0,XC1,YC0,YC1)
10.
KF1=withoutZero(Py*sp.diff(Phi1,y1),minDist)(XC0,XC1,YC0,YC1)
11.
MtrA=KF0*Nrml[0]+KF1*Nrml[1]
12.
k=0
13.
for i in range(len(MtrA)):
14.
for j in range(len(MtrA)):
15.
if(k>len(Akoef)-1):
16.
k=0
17.
MtrA[i,j]=(-2*lmd*MtrA[i,j]*h/2)*Akoef[k]
18.
k=k+1
19.
MtrA=MtrA+ ( np.eye(nrw) if isDiag else 0 )
20.
vtrf=2.*lmd*vphFun(q,x,z,Phi1)-2*alph*x[1]
21.
return MtrA,vtrf
Листинг 20 – Программная реализация модуля, который решает
интегральное (метод квадратурной формулы Гаусса) и дифференциальное
уравнение.
Строка 13 – 18 идет применение квадратурной формулы Гаусса. 19 –
находим левую часть уравнения, 20 – правую.
46
Рисунок 13 – График решения задачи математического моделирования
об опускании бугра грунтовых вод под действием силы тяжести при
800 точках, где интегральное уравнение решается квадратурной
формулой Гаусса.
На (рис. 13) мы видим, как с течением времени смещается кривая L. Но
по сложно сделать какие-то конкретные выводы, поэтому необходимо
заполнить таблицу получаемыми данными в ходе варьирования параметров
для этой задачи.
2.4.3. Исследование влияния численных методов
Проведем сопоставление результатов численного решения[13,14]
задачи
об
оседании
бугра
при
решении
интегрального
уравнения
47
квадратурной формулой Гаусса и формулой прямоугольника (рис 12,13).
Результаты представлены в табл. 1, 2. В данных таблице погрешность
|
, где T — время опускания бугра, которое мы варьируем, от
|
первоначальной высоты, полученное в ходе численного решения задачи.
предыдущее значение.
n
200
304
400
504
600
704
800
0,1
3,6
3,5
3,5
3,5
3,6
3,8
3,9
2,7
0
0
2,7
0,5
0,26
3,5
3,5
3,5
3,45
3,45
3,4
2,7
0
0
1,5
0
1,5
3,5
3,5
3,5
3,475
3,45
3,425
2,7
0
0
0,71
0,71
0,72
3,5
3,5
3,5
3,4625
3,4375
3,4125
3,1
0
0
1,1
0,73
0,72
3,125
3,125
3,125
3,125
3,125
3,125
0
0
0
0
0
0
,%
0.1/2.0
3,6
,%
0.1/4.0
3,6
,%
0.1/8.0
3,6125
,%
0.1/16.0
3,125
,% 0
Таблица 1. Влияние n и dt на погрешность при применении метода
квадратурной формулы Гаусса.
-
48
n
200
304
400
504
600
704
800
0,1
3,7
3,7
3,6
3,6
3,6
3,6
3,6
0
2,7
0
0
0
0
3,7
3,65
3,65
3,65
3,65
3,6
1,3
1,3
0
0
0
1,3
3,7
3,65
3,65
3,625
3,625
3,625
1,3
1,3
0
0,68
0
0
3,7
3,6625
3,65
3,6375
3.6375
3,625
1
1
0,3
0,3
0
0,3
3,7
3,65625 3,65
3.6375
3.6375
3,63125
1
1
0,34
0
0,17
,%
0.1/2.0
3,75
,%
0.1/4.0
3,75
,%
0.1/8.0
3,7375
,%
0.1/16.0
3,7375
,%
0,17
Таблица 2. Влияние n и dt на погрешность при применении формулы
прямоугольника.
В
строках
таблицы
указывается
время
опускания
бугра
до
определенного момента заданного программно, в нашем случае до середины,
в зависимости от количества узлов и dt(прирощения).
Анализируя таблицы, видим, что с уменьшением шага по времени и
длины отрезков, образующих подвижную границу, погрешность численного
решения уменьшается, т.е. численное решение сходится к точному.
Так же сравнивая таблицу 1 и 2, можно сделать вывод, что метод
49
Гаусса быстрее приближает решение задачи к его точному аналогу.
50
Заключение
В данной работе проводится исследование влияния применения
различных численных методов для решения интегрального уравнения на
примере задачи математического моделирования процесса оседания бугра
грунтовых вод в однородных слоях под действием силы тяжести, и
построения программного обеспечения[13,14] для решения проблемы
оседания бугра грунтовых водв однородных слоях с применением более
точной формулы для решения интегрального уравнения. Подобные задачи
являются востребованными, так как помогают преодолевать серьезные
экологические последствия затопления территорий, подтопления зданий и
сооружений, а применение более точных методов дает нам возможность
уменьшить погрешности решения, которое мы находим.
Очевидно, что для решения подобных задач численными методами
должны применять алгоритмы, обеспечивающие высокую степень точности
результата. Поэтому в данной работе я провожу исследования влияния
различных методов решения интегрального уравнения и делаю заключение,
какой из методов следует использовать при решении различных задач. Также
одним из современных методов приближенного решения является метод
дискретных особенностей, который дает возможность численно решить
задачу с высокой степенью достоверности результата. В данной работе
применяется метод дискретных особенностей для численного решения
дифференциальных уравнений.
По проведенным мной исследованиям, т.е решения интеграла,
интегрального уравнения с заданным точным решением вне задачи,
интегрального
уравнения,
которое
я
разбиваю
на
подотрезки
для
последующего решения более серьезных задач(где 8-16 отрезков не
достаточно для нахождения оптимального решения), и далее я использую
методы для решения задачи
математического моделирования процесса
оседания бугра грунтовых вод в однородных слоях под действием силы
51
тяжести, составляю таблицы по полученным мной данным, прихожу к
выводу, что самое близкое к точному решению дает квадратурная формула
Гаусса.
Достигнута следующая цель: разработано программное обеспечение
для решения интегральных уравнений квадратурной формулой Гаусса и
добавление его в программное обеспечение, которое реализует решение
задачи математического моделирования процесса оседания бугра грунтовых
вод в однородных слоях под действием силы тяжести.
Были выполнены следующие задачи:
1.
Изучена литература по рассматриваемому вопросу и
сделано заключение о необходимости и эффективности решения задачи
одним из численных методов;
2.
Для разработанной модели была получена численная схема
для реализации задачи различными методами;
3.
Разработан программный комплекс для численного счета на
одном из современных языков программирования;
4.
Исследовано полученное приближенное решение задачи и
проанализированы графики и таблицы, демонстрирующие поведение
различных методов решения интегрального уравнения в данной
задачи;
5.
Работа
Сделаны выводы по результатам исследований.
выполнена
на
языке
Python,
который
является
высокоуровневым языком программирования общего назначения. Так же
используются вспомогательные библиотеки языка, такие как NumPy и SimPy
[3] [4].
Работа имеет теоретическое и практическое значение, т.к. полученные
теоретические результаты могут быть использованы для корректировки
модели и в дальнейшем выбора метода дающего наиболее точное решение, а
разработанный комплекс программ может применяться для дальнейших
52
исследований и решений задач прикладного характера,
необходимо решать интегральные уравнения
в которых
53
Список литературы
1.
Никольский Д.Н., Дорофеева В.И. Математическое моделирование
двумерного процесса изменения уровня грунтовых вод под действием силы
тяжести методом дискретных особенностей// Вычислительные методы и
программирование,2011,12,c.85-89,http:// num-meth.srcc.msu.ru.
2.
Полубаринова-Кочина П.Я. Теория движения грунтовых вод. М.: Наука,
1977.
3.
Python 2.7.3 /https://www.python.org/ (Дата обращения: 12.04.2017)
4.
NumPy 1.8.1 /http://www.scipy.org/ (Дата обращения: 15.04.2017)
5.
Н. В. Копченова, И. А. Марон. Вычислительная математика в примерах
и задачах.-М.:Наука, 1972.
6.
Gnuplot 4.6 /http://www.gnuplot.info/ (Дата обращения: 6.05.2017)
7.
Голубева О.В. Курс механики сплошных сред.М.: Высшая школа, 1972.
8.
Чарный И.А. Подземная гидродинамика. М.: Гостоптехиздат, 1963. 396
с.
9.
Н. В. Копченова, И. А. Марон. Вычислительная математика в примерах
и задачах.-М.:Наука, 1972.
10.
Ракитин
В.А.,
методамвычислений
Первушин
с
В.Е.
приложением
Практическое
программ
руководство
для
по
персональных
компьютеров. М.:Высш. школа, 1998. 384 с.
11.
Макр Лутц, Программирование на Python, 4-е изд.–Символ-Плюс,2011
12.
Matplotlib 1.3.1 /http://www.matplotlib.org/ (Дата обращения: 23.04.2017)
13.
Голубев Г.В., Тумашев Г.Г. Фильтрация несжимаемой жидкости в
неоднородной пористой среде. Казань: Издательство КГУ, 1972. 196 с.
14. Пивень В.Ф. Математическое моделирование граничных задач гидродинамики в неоднородных слоях// Дис. ... докт. физ.-мат. наук: 05.13.18. Орел,
1998. 266 с.
15. Аксюхин А.А. Математическое моделирование граничных задач фильтрации к скважине в неоднородных слоях грунта. Дис. ... канд. физ.-мат.
54
наук: 05.13.18. Орел, 2000. 153с.
16. Никольский Д.Н. Математическое моделирование работы системы
скважин в однородных и неоднородных слоях с подвижной границей раздела
жидкостей различной вязкости. Дис. канд. физ.-мат. наук: 05.13.18. Орел,
2001. 191 с.
17. Ставцев С.Л. Исследование трехмерных граничных задач о дебите
системы несовершенных скважин в кусочно-неоднородных слоях. Дис....
канд. физ.- мат. наук: 05.13.01. М., 2003. 173 с.
18. Пивень В.Ф. Теоретическая механика, Орел, 2001.
19. Бухгольц Н.Н. Основной курс теоретической механики. М.:Наука,
1969,4.1,2.
20. Радыгин В.М., Голубева О.В. Применение функций комплексного переменного в задачах физики и техники. М.Высш. шк., 1983. 160 с.
55
Приложения
Листинг 21
Программная реализация метода прямоугольников для решения задачи
математического моделирования опускания бугра грунтовых вод под
действием силы тяжести.
1. #!/usr/bin/python
2. import numpy as np
3. import sympy as sp
4. from pylab import plot, show, grid, quiver,axis
5. from sys import path
6. path.append("./../")
7. from sfunctp import *
8.
9. from scipy.interpolate import interp1d
10.
11.# K=y**s, s=2; H=1
12.Kx=1
13.Hx=1
14.Py=1
15.d2 = 0.0
16.Phi1=-0.5/sp.pi*sp.log(sp.sqrt((x0-y0)**2+(x1-y1)**2))
17.#+0.5/sp.pi*sp.log(sp.sqrt((x0-(2*d2-y0))**2+(x1-(2*d2-y1))**2))
18.Psi2=0.5/sp.pi*sp.log(sp.sqrt((x0-y0)**2+(x1-y1)**2))
19.#+0.5/sp.pi*sp.log(sp.sqrt((x0-(2*d2-y0))**2+(x1-(2*d2-y1))**2))
20.sp.pprint(Phi1)
21.sp.pprint(Psi2)
22.n=800
23.#q=-np.array([np.pi])
24.q=-np.array([0.0])
25.z=np.array([[0.0],[1.0]],dtype='d')
26.lambdat=1.0
27.alphat=-1.0
28.dt=0.1/8.0
29.# ?
30.rSource=0.05
31.t_K=0.2
32.H_K=1.0
56
33.R_K=1.0
34.c_K=1.0
35.H0=H_K/np.pi*(np.arctan(R_K/(c_K*t_K))-np.arctan(-R_K/(c_K*t_K)))
36.
37.r_skv = H0/5.0
38.
39.H_min = H0/5.0+d2
40.print (H_min)
41.
42.thet=np.linspace(-20,20,n+1,endpoint=True)
43.x=np.array([thet,H_K/np.pi*(np.arctan((thet+R_K)/(c_K*t_K))np.arctan((thet-R_K)/(c_K*t_K)))+d2])
44.# x=np.array([tau+0.5*tau/(1+tau*tau),0.5/(1+tau*tau)+d2])
45.
46.plot(x[0],x[1])
47.#axis('equal')
48.#grid(True)
49.#show()
50.
51.
52.isMove=True
53.T=0
54.cycle=0
55.
56.while isMove:
57. A,f=ieFredgolm2dk(x,x,lambdat,alphat,q,z,varphiSources,Py,Phi1)
58. g=np.linalg.solve(A,f)
59.
60. W=VSources(q,x,z,Kx,Phi1)+VDL(x,x,g,Hx,Psi2)
61. #if cycle==0: quiver(x[0],x[1],W[0],W[1])
62.x+=W*dt
63.
64. x[0][0]=-20; x[0][n]=20
65. interpolant=interp1d(x[0].ravel(), x[1].ravel())
66.x1new=interpolant(thet)
67. x=np.array([thet,x1new])
68.
69.T+=dt; cycle+=1;
70. if cycle%50==0: print (cycle) ; plot(x[0],x[1])
71.
72.ymax=np.max(x[1])
73. print ("%12.10f" % T, ymax)
57
74. eps_H=0.00000000001
75. ymaxH=np.abs(ymax-H_min)
76.
77. #if cycle==1500: isMove=False; plot(x[0],x[1])
78.
79. if (ymaxH<=eps_H or ymax<=H_min ): isMove=False; plot(x[0],x[1])
80.
81.# print "%12.10f" % 1.4
82.
83.
84.print (cycle,"%12.10f" % T)
85.axis('equal')
86.grid(True)
87.show()
58
Листинг 22
Модуль реализующий метод прямоугольников.
1. #!/usr/bin/python
2.
3.
4.
5.
6.
7.
8.
9.
# agreements:
# border has the format plus point
# global symbol's variables:
#
x=(x0,x1) --- point collocation ( row )
#
y=(y0,y1) --- point with the singularity ( column )
# z --- point z contains a source (sink)
#
# a,b,c,d --- vector; A,B,C,D --- matrix
10.import numpy as np
11.import sympy as sp
12.myeps=0.0005
13.myeps2=myeps**2
14.x0,x1,y0,y1=sp.symbols('x0 x1 y0 y1')
15._Phi1_=-0.25/sp.pi*sp.log((x0-y0)**2+(x1-y1)**2)
16._Psi2_=-_Phi1_
17._Kx_=1
18._Hx_=1
19._Px_=1
20._Ky_=1
21._Hy_=1
22._Py_=1
23.def withoutZero(sFunct=None,minDist=myeps):
24.def vectorFunct(A,B,C,D):
i. gd=((A-C)**2+(B-D)**2)>minDist**2
ii. AA=A.copy()
iii. AA[np.logical_not(gd)]+=1
iv. result=np.zeros_like(A)
25.result[gd]=sp.lambdify((x0,x1,y0,y1),sFunct,modules='numpy')(AA,B,C,D
)[gd]
i. return result
26.return vectorFunct
59
27.def lmbd(sFunct=None):
28.def vectorFunct(A,B,C,D):
i. return
sp.lambdify((x0,x1,y0,y1),sFunct,modules='numpy')(A,B,C,D)
29.return vectorFunct
30.def varphiSources(q,x,z,Phi1=_Phi1_):
31.Z0,X0=np.meshgrid(z[0].ravel(),x[0].ravel())
32.Z1,X1=np.meshgrid(z[1].ravel(),x[1].ravel())
33.return -np.sum(q*lmbd(Phi1)(X0,X1,Z0,Z1),axis=1)
34.def VSources(q,x,z,Kx=_Kx_,Phi1=_Phi1_):
35.Z0,X0=np.meshgrid(z[0].ravel(),x[0].ravel())
36.Z1,X1=np.meshgrid(z[1].ravel(),x[1].ravel())
37.V0=lmbd(Kx*sp.diff(Phi1,x0))(X0,X1,Z0,Z1)
38.V1=lmbd(Kx*sp.diff(Phi1,x1))(X0,X1,Z0,Z1)
39.return -np.array([np.sum(q*V0,axis=1), np.sum(q*V1,axis=1)])
40.def Normal(x):
41.n=x[0].shape[0]-1
42.return np.array([-x[1,1:]+x[1,:n],x[0,1:]-x[0,:n]])
43.def normal(x):
44.nr=Normal(x)
45.return nr/np.sqrt(nr[0]**2+nr[1]**2)
46.def
ieFredgolm2dk(xrw,xcl,lmd,alph,q,z,vphFun,Py=_Py_,Phi1=_Phi1_,isDiag
=1,minDist=myeps):
47.nrw=xrw[0].shape[0]-1
48.ncl=xcl[0].shape[0]-1
49.x=np.array([(xrw[0,:nrw]+xrw[0,1:])/2,(xrw[1,:nrw]+xrw[1,1:])/2])
50.y=np.array([(xcl[0,:ncl]+xcl[0,1:])/2,(xcl[1,:ncl]+xcl[1,1:])/2])
51.YC0,XC0=np.meshgrid(y[0,:],x[0,:])
52.YC1,XC1=np.meshgrid(y[1,:],x[1,:])
53.Nrml=Normal(xcl)
54.KF0=withoutZero(Py*sp.diff(Phi1,y0),minDist)(XC0,XC1,YC0,YC1)
55.KF1=withoutZero(Py*sp.diff(Phi1,y1),minDist)(XC0,XC1,YC0,YC1)
56.MtrA=KF0*Nrml[0]+KF1*Nrml[1]
57.MtrA=-2*lmd*MtrA + ( np.eye(nrw) if isDiag else 0 )
58.vtrf=2.*lmd*vphFun(q,x,z,Phi1)-2*alph*x[1]
59.return MtrA,vtrf
60
60.def Theta(X0,X1,Y0,Y1,rEps=myeps):
61.R2=(X0-Y0)**2+(X1-Y1)**2
62.THETA=np.ones_like(R2)
63.rEps2=rEps**2
64.R2=R2/rEps2; R=R2**.5; R3=R2*R; R5=R3*R2; R7=R5*R2; R9=R7*R2
65.THETA[R2<rEps]=.125*(63.*R5-90.*R7+35.*R9)[R2<rEps]
66.return THETA
67.def V2(X0,X1,Y0,Y1,Hx=_Hx_,Psi2=_Psi2_,rEps=myeps):
68.return
Theta(X0,X1,Y0,Y1,rEps)*np.array([withoutZero(sp.diff(Psi2,x1)/Hx)(X0,
X1,Y0,Y1),withoutZero(-sp.diff(Psi2,x0)/Hx)(X0,X1,Y0,Y1)])
69.def MatrixVDL(xrw,xcl,Hx=_Hx_,Psi2=_Psi2_,rEps=myeps):
70.ncl=xcl[0].shape[0]-1
71.Y0,X0=np.meshgrid(xcl[0,:ncl],xrw[0].ravel())
72.Y1,X1=np.meshgrid(xcl[1,:ncl],xrw[1].ravel())
73.YNEXT0,X0=np.meshgrid(xcl[0,1:],xrw[0].ravel())
74.YNEXT1,X1=np.meshgrid(xcl[1,1:],xrw[1].ravel())
75.return
V2(X0,X1,Y0,Y1,Hx,Psi2,rEps)V2(X0,X1,YNEXT0,YNEXT1,Hx,Psi2,rEps)
76.def VDL(xrw,xcl,gg,Hx=_Hx_,Psi2=_Psi2_,rEps=myeps):
77.W=MatrixVDL(xrw,xcl,Hx,Psi2,rEps)*gg
78.return np.array([np.sum(W[0],axis=1),np.sum(W[1],axis=1)])
61
Листинг 23
Программная реализация квадратурной формулы Гаусса для решения
задачи математического моделирования опускания бугра грунтовых вод под
действием силы тяжести.
1.
2.
3.
4.
5.
6.
7.
8.
#!/usr/bin/python
import numpy as np
import sympy as sp
from pylab import plot, show, grid, quiver,axis
from sys import path
path.append("./../")
from prob import *
from scipy.interpolate import interp1d
9. Kx=1; Hx=1; Py=1
10.d2 = 0.0
11.Phi1=-0.5/sp.pi*sp.log(sp.sqrt((x0-y0)**2+(x1-y1)**2))
12.Psi2=0.5/sp.pi*sp.log(sp.sqrt((x0-y0)**2+(x1-y1)**2))
13.sp.pprint(Phi1)
14.sp.pprint(Psi2)
15.n=100
16.q=-np.array([0.0])
17.z=np.array([[0.0],[1.0]],dtype='d')
18.lambdat=1.0
19.alphat=-1.0
20.dt=0.1/8.0
21.# ?
22.rSource=0.05
23.t_K=0.2; H_K=1.0;
24.R_K=1.0; c_K=1.0
25.# H0=0.5
26.# r_skv = H0/5.0
27.# H_min = H0/5.0+d2
28.# print (H_min)
29.m = 8
30.mas_T = [-0.96028986,
-0.79666648,
-0.52553242,
-0.18343464,
62
0.18343464, 0.52553242, 0.79666648, 0.96028986]
31.mas_A = [0.1012854, 0.22238104, 0.31370664, 0.36268378, 0.36268378,
0.31370664, 0.22238104,0.1012854]
32.a = -20.; b = 20.
33.h = (b - a)/n; h1 = h/2.
34.x2 = a + h1;
35.k = 0; end=a
36.thet = np.arange(n*m)*1.
37.Akoef = np.arange(n*m)*1.
38.for i in range(n):
39.end=end+h
40.for j in range(m):
a. thet[k] = x2 + mas_T[j]*h1 + (i*h)
b. Akoef[k]=(end-(end-h))*(mas_A[j]/2)
c. k = k + 1
41.# print (thet)
42.x=np.array([thet,H_K/np.pi*(np.arctan((thet+R_K)/(c_K*t_K))np.arctan((thet-R_K)/(c_K*t_K)))+d2])
43.Akoef = Akoef[:len(Akoef)-1]
44.H0=H_K/np.pi*(np.arctan((0.+R_K)/(c_K*t_K))-np.arctan((0.R_K)/(c_K*t_K)))+d2
45.r_skv = H0/5.0
46.H_min = H0/5.0+d2
47.plot(x[0],x[1])
48.isMove=True
49.T=0; cycle=0
50.while isMove:
51.A,f=ieFredgolm2dk(x,x,Akoef,h1,lambdat,alphat,q,z,varphiSources,Py,Phi1
)
52.g=np.linalg.solve(A,f)
53.W=VSources(q,x,z,Kx,Phi1)+VDL(x,x,g,Hx,Psi2)
54.# if cycle==0: quiver(x[0],x[1],W[0],W[1])
55.x_new=x + W*dt
56.x_new[0][0]=x[0][0]; x_new[0][n*m-1]=x[0][n*m-1]
57.interpolant=interp1d(x_new[0].ravel(), x_new[1].ravel())
58.x1new=interpolant(x[0])
63
59.x=np.array([x[0],x1new])
60.T+=dt; cycle+=1;
61.if cycle%50==0: print (cycle) ; plot(x[0],x[1])
62.ymax=np.max(x[1])
63.print ("%12.10f" % T, ymax)
64.eps_H=0.00000000001
65.ymaxH=np.abs(ymax-H_min)
66.if cycle==500: isMove=False; plot(x[0],x[1])
67.if (ymaxH<=eps_H or ymax<=H_min ): isMove=False; plot(x[0],x[1])
68.print (cycle,"%12.10f" % T)
69.axis('equal')
70.grid(True)
71.show()
64
Листинг 24
Модуль реализующий метод квадратурной формулы Гаусса.
1. #!/usr/bin/python
2.
3.
4.
5.
6.
7.
8.
9.
# agreements:
# border has the format plus point
# global symbol's variables:
#
x=(x0,x1) --- point collocation ( row )
#
y=(y0,y1) --- point with the singularity ( column )
# z --- point z contains a source (sink)
#
# a,b,c,d --- vector; A,B,C,D --- matrix
10.import numpy as np
11.import sympy as sp
12.myeps=0.0000005
13.myeps2=myeps**2
14.x0,x1,y0,y1=sp.symbols('x0 x1 y0 y1')
15._Phi1_=-0.25/sp.pi*sp.log((x0-y0)**2+(x1-y1)**2)
16._Psi2_=-_Phi1_
17._Kx_=1
18._Hx_=1
19._Px_=1
20._Ky_=1
21._Hy_=1
22._Py_=1
23.def withoutZero(sFunct=None,minDist=myeps):
24.def vectorFunct(A,B,C,D):
i. gd=((A-C)**2+(B-D)**2)>minDist**2
ii. AA=A.copy()
iii. AA[np.logical_not(gd)]+=1
iv. result=np.zeros_like(A)
25.result[gd]=sp.lambdify((x0,x1,y0,y1),sFunct,modules='numpy')(AA,B,C,D)
[gd]
i. return result
26.return vectorFunct
27.def lmbd(sFunct=None):
65
28.def vectorFunct(A,B,C,D):
i. return
sp.lambdify((x0,x1,y0,y1),sFunct,modules='numpy')(A,B,C,D)
29.return vectorFunct
30.def varphiSources(q,x,z,Phi1=_Phi1_):
31.Z0,X0=np.meshgrid(z[0].ravel(),x[0].ravel())
32.Z1,X1=np.meshgrid(z[1].ravel(),x[1].ravel())
33.return -np.sum(q*lmbd(Phi1)(X0,X1,Z0,Z1),axis=1)
34.def VSources(q,x,z,Kx=_Kx_,Phi1=_Phi1_):
35.Z0,X0=np.meshgrid(z[0].ravel(),x[0].ravel())
36.Z1,X1=np.meshgrid(z[1].ravel(),x[1].ravel())
37.V0=lmbd(Kx*sp.diff(Phi1,x0))(X0,X1,Z0,Z1)
38.V1=lmbd(Kx*sp.diff(Phi1,x1))(X0,X1,Z0,Z1)
39.return -np.array([np.sum(q*V0,axis=1), np.sum(q*V1,axis=1)])
40.def Normal(x):
41.n=x[0].shape[0]-1
42.return np.array([-x[1,1:]+x[1,:n],x[0,1:]-x[0,:n]])
43.def normal(x):
44.nr=Normal(x)
45.return nr/np.sqrt(nr[0]**2+nr[1]**2)
46.def
ieFredgolm2dk(xrw,xcl,Akoef,h,lmd,alph,q,z,vphFun,Py=_Py_,Phi1=_Phi1
_,isDiag=1,minDist=myeps):
47.nrw=xrw[0].shape[0]-1
48.ncl=xcl[0].shape[0]-1
49.x=np.array([(xrw[0,:nrw]),(xrw[1,:nrw])])
50.y=np.array([(xcl[0,:ncl]),(xcl[1,:ncl])])
51.YC0,XC0=np.meshgrid(y[0,:],x[0,:])
52.YC1,XC1=np.meshgrid(y[1,:],x[1,:])
53.Nrml=Normal(xcl)
54.KF0=withoutZero(Py*sp.diff(Phi1,y0),minDist)(XC0,XC1,YC0,YC1)
55.KF1=withoutZero(Py*sp.diff(Phi1,y1),minDist)(XC0,XC1,YC0,YC1)
56.MtrA=KF0*Nrml[0]+KF1*Nrml[1]
57.# MtrA=-2*lmd*Akoef*h*MtrA+ ( np.eye(nrw) if isDiag else 0 )
58.# vtrf=2.*lmd*vphFun(q,x,z,Phi1)-2*alph*x[1]
59.k=0
60.for i in range(len(MtrA)):
66
i. for j in range(len(MtrA)):
ii. if(k>len(Akoef)-1):
1. k=0
iii. MtrA[i,j]=(-2*lmd*MtrA[i,j]*h/2)*Akoef[k]
iv. k=k+1
61.MtrA=MtrA+ ( np.eye(nrw) if isDiag else 0 )
62.vtrf=2.*lmd*vphFun(q,x,z,Phi1)-2*alph*x[1]
63.return MtrA,vtrf
64.def Theta(X0,X1,Y0,Y1,rEps=myeps):
65.R2=(X0-Y0)**2+(X1-Y1)**2
66.THETA=np.ones_like(R2)
67.rEps2=rEps**2
68.R2=R2/rEps2; R=R2**.5; R3=R2*R; R5=R3*R2; R7=R5*R2; R9=R7*R2
69.THETA[R2<rEps]=.125*(63.*R5-90.*R7+35.*R9)[R2<rEps]
70.return THETA
71.def V2(X0,X1,Y0,Y1,Hx=_Hx_,Psi2=_Psi2_,rEps=myeps):
72.return
Theta(X0,X1,Y0,Y1,rEps)*np.array([withoutZero(sp.diff(Psi2,x1)/Hx)(X0,
X1,Y0,Y1),withoutZero(-sp.diff(Psi2,x0)/Hx)(X0,X1,Y0,Y1)])
73.def MatrixVDL(xrw,xcl,Hx=_Hx_,Psi2=_Psi2_,rEps=myeps):
74.ncl=xcl[0].shape[0]-1
75.Y0,X0=np.meshgrid(xcl[0,:ncl],xrw[0].ravel())
76.Y1,X1=np.meshgrid(xcl[1,:ncl],xrw[1].ravel())
77.YNEXT0,X0=np.meshgrid(xcl[0,1:],xrw[0].ravel())
78.YNEXT1,X1=np.meshgrid(xcl[1,1:],xrw[1].ravel())
79.return
V2(X0,X1,Y0,Y1,Hx,Psi2,rEps)V2(X0,X1,YNEXT0,YNEXT1,Hx,Psi2,rEps)
80.def VDL(xrw,xcl,gg,Hx=_Hx_,Psi2=_Psi2_,rEps=myeps):
81.W=MatrixVDL(xrw,xcl,Hx,Psi2,rEps)*gg
82.return np.array([np.sum(W[0],axis=1),np.sum(W[1],axis=1)])
1/--страниц
Пожаловаться на содержимое документа