Численное решение краевых задач математической физики методом сеток
Федеральное
государственное бюджетное
высшего профессионального образования
«САМАРСКИЙ ГОСУДАРСТВЕННЫЙ
УНИВЕРСИТЕТ имени академика С.П. КОРОЛЕВА
(национальный
Факультет информатики
Кафедра технической кибернетики
Расчётно-пояснительная записка к курсовому проекту
по дисциплине «Численные методы»
Тема: «ЧИСЛЕННОЕ РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ МАТЕМАТИЧЕСКОЙ ФИЗИКИ
МЕТОДОМ СЕТОК»
Вариант № 27
Выполнил студент
Группа 6409
Руководитель работы
2 0 1
2
Задание
Осуществить математическую постановку краевой задачи для физического процесса, описанного в предложенном варианте курсового проекта.
Осуществить построение сеточной схемы, приближающей полученную краевую задачу. При этом следует согласовать с преподавателем тип сеточной схемы.
Провести аналитическое исследование схемы: показать, что схема аппроксимирует исходную краевую задачу, и найти порядки аппроксимации относительно шагов дискретизации; исследовать устойчивость схемы и сходимость сеточного решения к решению исходной задачи математической физики.
Разработать (описать) алгоритм численного решения сеточной краевой задачи.
Разработать компьютерную программу,
реализующую созданный
Используя разработанную
компьютерную программу провести исследование
погрешности сеточного решения
методом вычислительного
Оформить пояснительную записку к курсовой работе в соответствии с требованиями, предъявляемыми к учебным текстовым документам (СТО СГАУ 02068410-004-2007, Общие требования к учебным текстовым документам. – Самара, СГАУ, 2007).
Вариант 27
В тонком однородном стержне длиной и сечением происходит тепловой процесс на временном промежутке . Боковая поверхность стержня теплоизолированная, теплоизолированными также являются и концы стержня. В начальный момент времени температура стержня описывалась функцией . Коэффициенты теплопроводности и объемной теплоемкости материала стержня равны и соответственно.
Для численного решения использовать конечно-разностную схему Кранка-Николсона.
Значения параметров:
Реферат
Пояснительная записка: 25 с., 7 рисунков, 4 источника.
УРАВНЕНИЯ МАТЕМАТИЧЕСКОЙ ФИЗИКИ, ТЕПЛОВОЙ ПРОЦЕСС, АППРОКСИМАЦИЯ КРАЕВОЙ ЗАДАЧИ КОНЕЧНО-РАЗНОСТНОЙ СХЕМОЙ, РАЗНОСТНАЯ СХЕМА КРАНКА-НИКОЛСОНА, АППРОКСИМАЦИЯ, УСТОЙЧИВОСТЬ, СХОДИМОСТЬ, СКОРОСТЬ СХОДИМОСТИ.
Основным объектом исследования
данного курсового проекта
Целью данной работы
является создание компьютерной программы
для расчета изменения
Для решения поставленной задачи строится разностная схема Кранка-Николсона. Для схемы найден порядок аппроксимации, показана условная устойчивость и выполнение необходимого признака Неймана. Для реализации разностной схемы спроектирован программный модуль, позволяющий получить результаты заданными параметрами построения сетки. Проведена работа по исследованию качества полученного решения с помощью вычислительного эксперимента: сопоставление аналитического и сеточного решений, исследование скорости сходимости.
Программа написана на Java: 1.6.0_26, использовалась виртуальная машина Java HotSpot(TM) Client VM 20.1-b02.
Содержание
Задание 2
Реферат 4
Содержание 5
Введение 6
1 Постановка краевой задачи 7
2 Построение разностной схемы Кранка-Николсона 8
3 Аппроксимация 11
4 Устойчивость 13
5 Необходимый признак Неймана 15
6 Моделирование процесса на компьютере с помощью разностной схемы 16
7 Экспериментальное исследование скорости сходимости 19
Заключение 24
Список использованных источников 25
Приложение 26
Введение
Спектр краевых задач математической физики, для которых возможно получить аналитическое решение, довольно узкий. Одним из подходов, в рамках которого можно получить решение гораздо более широкого множества краевых задач математической физики, является метод конечных разностей. Этот метод позволяет осуществить построение разностных схем, в данной работе будет использована схема Кранка-Николсона.
Для решения краевой задачи сначала строится разностная схема методом замены производных разностными отношениями. Далее проводится исследование схемы, а именно исследование порядков аппроксимации численного решения на регулярной сетке, а также исследование устойчивости решения, удовлетворения признаку устойчивости Неймана.
На основе полученной разностной
схемы разработана программа, и
в данной работе будут приведены
некоторые полученные результаты в
сравнении с аналитическим
1 Постановка краевой задачи
Уравнение теплопроводности [1] для стержня при отсутствии бокового теплообмена имеет вид .
Отсутствие теплообмена на концах стержня выражается в виде условий:
Добавляя к этим краевым условиям начальное условие , получим задачу математической физики:
(1)
Введем обозначения:
– линейный оператор, действующий на функцию . Результат действия оператора на пусть имеет вид:
.
Также рассмотрим вектор правых частей уравнений:
.
Если соблюдены ограничения на переменные из системы (1), то краевая задача (1) представима в виде:
,
2 Построение разностной схемы Кранка-Николсона
Определим равномерную сетку как множество узлов :
, где
Сеточная область является аналогом области , для которой в разделе 1 поставлена краевая задача. На сеточной области определены сеточные функции, образующие пространство сеточных функций .
Построение схемы Кранка-
Рисунок 1 – Шаблон схемы Кранка-Николсона
(
Обозначим значение функции в узле :
.
Рассмотрим наряду с основными узлами шаблона вспомогательный узел , обозначенный крестом на рисунке 1. Значение функции нем обозначим:
.
Выразим производные, входящие в состав системы (1), через разностные производные, получим приближенные равенства:
, (2.1)
, (2.2)
. (2.3)
Запишем первое уравнение системы (1) для узла и :
, (2.4)
(2.5).
Далее подставим (2.1) и (2.2) в (2.4); также подставим (2.1) и (2.3) в (2.5), увеличив в уравнении (2.1) индекс на единицу. Получим следующую систему:
(2.6)
где - значение решения во вспомогательном узле.
Для производной на границах используем разностную аппроксимацию производной первого порядка (соответственно левую для правой границы и правую для левой границы):
Теперь если умножить каждое из уравнений (2.6) на число 0,5 и сложить, а также дописать начальное и граничные условия, то получим схему Кранка-Николсона для нашей задачи:
(2.7)
где .
Эта схема является неявной и для расчёта решения мы используем метод прогонки.
Введем обозначения:
– линейный оператор, действующий на сеточную функцию , , значения которой в узлах с соответствующими номерами равны . Результат действия оператора на пусть имеет вид:
.
Также рассмотрим вектор – одноиндексную сеточную функцию:
, .
Тогда схема (2.7) представима в виде:
.
3 Аппроксимация
Рассмотрим общий вид краевой задачи относительно функции :
Здесь - линейный оператор, а – составной функциональный объект.
Для решения этой задачи мы используем ее сеточный аналог:
Сеточное решение
где
Сначала рассмотрим Разложим каждое слагаемое в выражении невязки в окрестностях некоторой точки Для простоты введем обозначение Также обозначим как Будем действовать поэтапно. Выпишем отдельно
Теперь найдем разложение в ряд первого слагаемого:
Во втором слагаемом нам понадобится
Теперь мы можем найти разложение выражения
Получим разложение числителя второго слагаемого:
И выражение числителя первого слагаемого:
Окончательно, используя ранее полученные промежуточные результаты, узнаем порядок локальной аппроксимации первого сеточного уравнения:
Как видно из локальной аппроксимации первого разностного уравнения в точке всилу произвольности выбора индексов
Проделаем аналогичный путь для :
Как видим, второе разностное уравнение аппроксимирует второе уравнение задачи (1) со сколь угодно высоким порядком.
Рассмотрим далее :
Локальный порядок аппроксимации:
Аналогично :
Локальный порядок аппроксимации:
В целом для всей системы уравнений можно сказать, что схема аппроксимирует исходную задачу, и погрешность аппроксимации схемы .
4 Устойчивость
Чтобы разностная схема была устойчивой, в первую очередь необходимо существование решения и притом единственного, также необходимо выполнение следующего условия:
где - возмущение правой части (погрешность схемы, неточность модели, вычислительная погрешность) [3].
Учитывая то, что наша схема линейна, задачу нахождения отклика на возмущение можно записать следующим образом:
(4.1)
Определим нормы в пространствах и таким образом:
где и - возмущения правой части и .
Пусть , перепишем уравнение (4.1):
Далее предположим, что и воспользуемся неравенством треугольника:
Раз уж записанное выше равенство справедливо, то и подавно справедливо будет
Заметим, что правая часть неравенства не зависит от коэффициента , а значит, справедлива для всех значений коэффициента, и тогда, попутно приведя подобные слагаемые, запишем:
Нам интересно, ограничена ли норма относительно возмущения правой части, для того чтобы это узнать, рассмотрим последнее уравнение при различных значениях
…
Учтем, что , и получим
Итак, мы доказали, что схема устойчива при выполнении условия .
5 Необходимый признак Неймана
Необходимый признак Неймана для разностных задач определяет необходимое условие устойчивости схемы. Таким образом, если он не выполняется, то схема не является устойчивой [4].
Теорема звучит так: для устойчивости разностной схемы (2.7) необходимо, чтобы существовали такие постоянные и , что для любых сеток мелкостью и для любых значений параметра было справедливо неравенство .
Собственные числа оператора будем искать, сделав замену , и сделав в нашей схеме (без граничных условий) замену
где -мнимая единица.
Первое уравнение системы с учетом замены примет вид
Найдем выражение для собственного числа
Заметим, что , кроме того , значит
Что говорит нам о выполнении условия устойчивости Неймана.
6 Моделирование процесса на компьютере с помощью разностной схемы
Получено соотношение, которое позволяет реализовать разностную схему программно:
Из начальных условий нам известны значения функции на временном слое , а также соотношения на граничных слоях и Этого достаточно, чтобы определить значения функции во всех узлах на всех временных слоях.
Для получения значений функции на каждом временном слое будем решать систему линейных уравнений , где матрицы и вектора имеют следующий вид
где
Для решения будем использовать метод прогонки, обладающий асимптотической сложностью O(n). Код программы, использовавшейся для моделирования процесса, приведен в приложении.
Ниже приведем некоторые результаты на графиках, отображающих решение, полученное разностным методом и аналитическим методом решения.
Рисунок
2 – График решения при фиксированной
координате
Рисунок
3 – График решения при фиксированной
координате
Рисунок
4 – График решения при фиксированном
времени
Рисунок
5 – График решения при фиксированном
времени
7 Экспериментальное исследование скорости сходимости
Для экспериментального исследования скорости сходимости явной схемы изучим скорость изменения погрешности при уменьшении шагов сетки.
Зададим погрешность сеточного решения относительно аналитического следующим образом:
где и - некоторые константы, а .
Спланируем эксперимент следующим образом:
Зафиксируем определенное количество шагов сетки (по координате или по времени), а по другой оси его будем увеличивать. На каждой итерации процесса будем вычислять получившуюся погрешность. Погрешность будем рассчитывать как разницу с аналитическим решением в заданной точке, посчитанной с точностью
Результаты представим на графиках, чтобы визуально оценить характер сходимости, а также в таблицах.
Зафиксируем 400 разбиений по пространственному измерению, получим
Рисунок
6– Зависимость погрешности решения от
мелкости сетки по
Приведем таблицу значений погрешностей для этого случая
|
K |
I |
|
|
|
|
10 |
400 |
3.0 |
0.058 |
|
20 |
400 |
1.5 |
0.041 |
1.41 |
40 |
400 |
0.75 |
0.029 |
1.41 |
80 |
400 |
0.375 |
0.0205 |
1.41 |
160 |
400 |
0.188 |
0.0146 |
1.40 |
320 |
400 |
0.094 |
0.0104 |
1.36 |
640 |
400 |
0.047 |
0.0075 |
1.389 |
1280 |
400 |
0.023 |
0.00557 |
1.35 |
2560 |
400 |
0.012 |
0.0042 |
1.308 |
5120 |
400 |
0.006 |
0.00335 |
1.272 |
Таблица 1 – исследование скорости сходимости.
Для анализа зависимости от величины шага по пространственной координате, будем рассчитывать количество шагов по временному измерению по эвристической формуле чтобы гарантировать, что погрешность, получаемая за счет величины шага по времени, заведомо меньше, чем результирующая. Приведем график и таблицу результатов.
Рисунок
7– Зависимость погрешности решения от
мелкости сетки по
|
K |
I |
|
|
|
|
600 |
10 |
0.3 |
0.249 |
|
2100 |
20 |
0.15 |
0.125 |
1.992 |
8100 |
40 |
0.075 |
0.062 |
2.01 |
32100 |
80 |
0.038 |
0.028 |
2.214 |
Таблица 2 – исследование скорости сходимости.
По итогам эксперимента можно сделать вывод о справедливости нашей аналитической оценки скорости сходимости от величины шага по измерению координаты, но с зависимостью от величины шага по времени не все очевидно. По данным в таблице видно, что скорость сходимости не достигает полученных аналитически оценок, и, мало того, убывает при измельчении сетки.
Заключение
В процессе выполнения работы методом конечных разностей была построена разностная схема Кранка-Николсона для исходной задачи математической физики. Для схемы было проведено исследование аппроксимации, была установлена условная устойчивость и выполнение необходимого признака Неймана.
Был создан программный модуль, моделирующий процесс распространения тепла в тонком стержне. Для проверки правильности подсчета порядка аппроксимации был проведён вычислительный эксперимент с целью установить зависимость погрешности численного решения задачи от шага сетки. При моделировании была подтверждена сходимость решения к истинному на рассмотренных сетках, однако остался незакрытым вопрос о скорости этой сходимости при измельчении шага по времени, т.к. при эксперименте аналитическая оценка не подтвердилась.
В результате были выделены следующие достоинства метода конечных разностей:
- широкая распространенность (большое количество доступных источников информации);
- сходимость к аналитическому решению (в данной работе она наблюдалась);
- подходит для решения класса задач, исследованного в данной работе;
- простота построения дискретных моделей.
Отметим, что конечных разностей не лишён следующих недостатков:
- сложность теоретического исследования построенной модели и доказательства её сходимости к исходной задаче;
- Довольно невысокая точность и большая вычислительная сложность по сравнению с аналитическими методами решения данного класса задач.
Список использованных источников
- А.А. Самарский Теория разностных схем. Учебное пособие. – М.: Наука, 1983. – 616 с.
- Е.Н. Сметанникова Конспект лекций по курсу «Численные методы математического анализа». – Самара: Самарский государственный аэрокосмический университет. 1999. – 97с.
- Дегтярев А.А. Метод конечных разностей. Электронное учебное пособие. – Самара: Самарский государственный аэрокосмический университет. 2011. – 83с.
- Дегтярев А.А. Численные методы решения задач математической физики. – Самара: Самарский государственный аэрокосмический университет. 2009. – 26с.
Приложение
22 public class GridComputer {
23
24 private static int big_i = 10000;
25 private static int big_k = 10000;
26 private static double k = 0.067;
27 private static double c = 1.84;
28 private static double length = 3;
29 private static double big_t = 30;
30 private static double x;
31 private static double t;
32
41 private static double maxDifference(double[][] numerical, int big_i, int big_k, double h_x, double t) {
42 double result = Double.NEGATIVE_INFINITY;
43 double current = 0;
44 SeriesComputer computer =
new SeriesComputer(length, k, c, 100500, big_i + 1, big_k + 1, 0, 0, big_t);
45 for (int i = 0; i < numerical[1].length; i++) {
46 current = Math.abs(numerical[1][i] - computer.computeInPoint(h_x * i, t, 30, big_t));
47 if (current > result){
48 result = current;
49 }
50 }
51 return result;
52 }
53
54 // Максимальное расхождение
сеточного и аналитического
55 private static double eps(int big_i, int big_k) {
56 double result = Double.NEGATIVE_INFINITY;
57 double current;
58
59
60 double[][] currentLayer = new double[2][big_i + 1];
61 double h_t = big_t / big_k;
62 double h_x = length / big_i;
63 for (int i = 0; i < big_i + 1; i++) {
64 currentLayer[0][i] = i * h_x;
65 currentLayer[1][i] = initialFunction(i * h_x);
66 }
67 double gamma = k / c * h_t / (h_x * h_x) * 0.5;
68 for (int k = 1; k <= big_k; k++) {
69
currentLayer[1] = computeOneLayer(currentLayer[
70 if (k > 1) {
71 current = maxDifference(currentLayer, big_i, big_k, h_x, h_t * k);
72
73 if (result < current) {
74 result = current;
75 }
76 }
77 }
78
79
80 return result;
81 }
82
180
181 // значение t должно принадлежать сетке
182 public static double[][] computeSeriesWhenT() {
183 double h_t = big_t / big_k;
184 long layer = Math.round(t / h_t);
185 double[][] result = computeLayer(layer);
186 return result;
187 }
188
189 public static double[][] computeSeriesWhenX() {
190 double h_t = big_t / big_k;
191 double h_x = length / big_i;
192 double[][] result = new double[2][big_k + 1];
193 int target_i = (int) Math.round(x / h_x);
194 double[][] currentLayer = new double[2][big_i + 1];
195 for (int i = 0; i < big_i + 1; i++) {
196 currentLayer[0][i] = i * h_x;
197 currentLayer[1][i] = initialFunction(i * h_x);
198 }
199 result[0][0] = 0;
200 result[1][0] = currentLayer[1][target_i];
201 double gamma = k / c * h_t / (h_x * h_x) * 0.5;
202 for (int k = 1; k <= big_k; k++) {
203
currentLayer[1] = computeOneLayer(currentLayer[
204 result[0][k] = k * h_t;
205 result[1][k] = currentLayer[1][target_i];
206 }
207 return result;
208 }
209
210 public static double initialFunction(double x) {
211 return Math.sin(Math.PI * x / length) + 10;
212 }
213
214 /**
215 * изображает трехдиагональную матрицу тремя массивами - диагоналями внутри
216 */
217 public static class ThreeDiagMatrix {
218
219 private double[] bot;
220 private double[] mid;
221 private double[] top;
222
223 public int size() {
224 return mid.length;
225 }
226
227 public double get(int i, int j) {
228 if (Math.abs(i - j) > 1) {
229 return 0;
230 } else if (i == j) {
231 return mid[j];
232 } else { // |i-j|=1
233 if (i - j == 1) {
234 return bot[j];
235 } else { // (i - j == -1)
236 return top[i];
237 }
238 }
239 }
240
241 public void set(int i, int j, double value) {
242 if (Math.abs(i - j) > 1) {
243 return;
244 } else if (i == j) {
245 mid[j] = value;
246 } else { // |i-j|=1
247 if (i - j == 1) {
248 bot[j] = value;
249 } else { // (i - j == -1)
250 top[i] = value;
251 }
252 }
253 }
254
255 public ThreeDiagMatrix(double[] bot, double[] mid, double[] top) {
256 if (!((bot.length == top.length) && (bot.length + 1 == mid.length))) {
257
throw new java.lang.
258 }
259 this.bot = Arrays.copyOf(bot, bot.length);
260 this.mid = Arrays.copyOf(mid, mid.length);
261 this.top = Arrays.copyOf(top, top.length);
262 }
263 }
264
265 /*
266 * решает систему с трехдиагональной матрицей и столбцом справа b
267 */
268 public static double[] linsolve(ThreeDiagMatrix matrix, double[] b) {
269 int N = matrix.size() - 1;
270 // прямой ход
271 // делим первую строку
272 matrix.set(0, 1,
273 matrix.get(0, 1) / matrix.get(0, 0));
274 b[0] = b[0] / matrix.get(0, 0);
275 matrix.set(0, 0, 1.0);
276 // делим остальные строки
277 for (int i = 1; i < matrix.size() - 1; i++) {
278 matrix.set(i, i + 1,
279 matrix.get(i, i + 1) / (matrix.get(i, i) - matrix.get(i, i - 1) * matrix.get(i - 1, i)));
280 b[i] = (b[i] - matrix.get(i, i - 1) * b[i - 1]) / (matrix.get(i, i) - matrix.get(i, i - 1) * matrix.get(i - 1, i));
281 matrix.set(i, i - 1, 0);
282 matrix.set(i, i, 1.0);
283 }
284 // делим последнюю строку
285 b[N] = (b[N] - matrix.get(N, N - 1) * b[N - 1]) / (matrix.get(N, N) - matrix.get(N, N - 1) * matrix.get(N - 1, N));
286 matrix.set(N, N, 1.0);
287 matrix.set(N, N - 1, 0);
288 // обратный ход
289 double[] x = new double[matrix.size()];
290 x[N] = b[N];
291 for (int i = N - 1; i >= 0; i--) {
292 x[i] = b[i] - matrix.get(i, i + 1) * x[i + 1];
293 }
294 return x;
295 }
296
297 /*
298 * находит один слой по времени с помощью предыдущего.
299 * gamma = a/2 * h_t / h_x^2 , где a это коэффициент при второй производной
300 */
301 public static double[] computeOneLayer(double[] preLayer, double gamma) {
302 int size = preLayer.length;
303 double[] top = new double[size - 3];
304 double[] mid = new double[size - 2];
305 double[] bot = new double[size - 3];
306 double[] b = new double[size-2];
307 for (int i = 0; i < size-3; i++) {
308 bot[i] = -1 * gamma;
309 mid[i] = (i==0) ? (1 + gamma) : (1 + 2 * gamma);
310 top[i] = -1 * gamma;
311 b[i] = (i == 0)
312 ? ((1 - gamma) * preLayer[1] + gamma * preLayer[2])
313 : ((1 - 2 * gamma) * preLayer[i+1] + gamma * (preLayer[i] + preLayer[i+2]));
314 }
315 mid[size - 3] = 1 + gamma;
316 b[size - 3] = ((1 - gamma) * preLayer[size-2] + gamma * preLayer[size-3]);
317 ThreeDiagMatrix matrix = new ThreeDiagMatrix(bot, mid, top);
318 double[] out = linsolve(matrix, b);
319 double[] result = new double[size];
320 for (int i=0; i<size-2; i++) {
321 result[i+1] = out[i];
322 }
323 result[0] = result[1];
324 result[size-1] = result[size-2];
325 return result;
326 }
327
328 public static double[][] computeLayer(long target_k) {
329 double[][] currentLayer = new double[2][big_i + 1];
330 double h_t = big_t / big_k;
331 double h_x = length / big_i;
332 for (int i = 0; i < big_i + 1; i++) {
333 currentLayer[0][i] = i * h_x;
334 currentLayer[1][i] = initialFunction(i * h_x);
335 }
336 double gamma = k / c * h_t / (h_x * h_x) * 0.5;
337 for (int k = 1; k <= target_k; k++) {
338
currentLayer[1] = computeOneLayer(currentLayer[
339 }
340 return currentLayer;
341 }
342
452 }
