12.8. Метод конечных разностей решения краевых задач для оду

Постановка задачи: найти решение линейного дифференциального уравнения

, (1)

удовлетворяющего краевым условиям:. (2)

Теорема. Пусть . Тогда существует единственное решение поставленной задачи.

К данной задаче сводится, например, задача об определении прогибов балки, которая на концах опирается шарнирно.

Основные этапы метода конечных разностей:

1) область непрерывного изменения аргумента ([a,b]) заменяется дискретным множеством точек, называемых узлами: .

2) Искомая функция непрерывного аргумента x, приближенно заменяется функцией дискретного аргумента на заданной сетке, т.е. . Функция называется сеточной.

3) Исходное дифференциальное уравнение заменяется разностным уравнением относительно сеточной функции. Такая замена называется разностной аппроксимацией.

Таким образом, решение дифференциального уравнения сводится к отысканию значений сеточной функции в узлах сетки, которые находятся из решения алгебраических уравнений.

Аппроксимация производных.

Для аппроксимации (замены) первой производной можно воспользоваться формулами:

- правая разностная производная,

- левая разностная производная,

- центральная разностная производная.

т.е., возможно множество способов аппроксимации производной.

Все эти определения следуют из понятия производной как предела: .

Опираясь на разностную аппроксимацию первой производной можно построить разностную аппроксимацию второй производной:

(3)

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

Определение. Погрешностью аппроксимации n- ой производной называется разность: .

Для определения порядка аппроксимации используется разложение в ряд Тейлора.

Рассмотрим правую разностную аппроксимацию первой производной:

Т.е. правая разностная производная имеет первый по h порядок аппроксимации.

Аналогично и для левой разностной производной.

Центральная разностная производная имеет второй порядок аппроксимации.

Аппроксимация второй производной по формуле (3) также имеет второй порядок аппроксимации.

Для того чтобы аппроксимировать дифференциальное уравнение необходимо в нем заменить все производные их аппроксимациями. Рассмотрим задачу (1), (2) и заменим в(1) производные:

.

В результате получим:

(4)

Порядок аппроксимации исходной задачи равен 2, т.к. вторая и первая производные заменены с порядком 2, а остальные – точно.

Итак, вместо дифференциальных уравнений (1), (2) получена система линейных уравнений для определения в узлах сетки.

Схему можно представить в виде:

т.е., получили систему линейных уравнений с матрицей:

Данная матрица является трехдиагональной, т.е. все элементы, которые расположены не на главной диагонали и двух прилегающих к ней диагоналях равны нулю.

Решая полученную систему уравнений, мы получим решение исходной задачи.

Для решения таких СЛАУ имеется экономичный метод прогонки.

Рассмотрим метод прогонки для СЛАУ:

(1)

Решение данной системы ищем в виде:

(2)

Подставляя в первое уравнение, получим:

Здесь учтено, что данное соотношение должно выполняться при любом

Так как

, (3)

то подставляя (3) во второе уравнение, получим:

Сравнивая с (2) получим

.

Таким образом, можно найти все .

Тогда из последнего уравнения (1) находим:

Затем последовательно находим:

Таким образом, алгоритм метода прогонки можно представить в виде:

1) Находим

2) Для i=1,n-1: (4)

3) Находим

4) Для i=n-1 до 1 находим:

Шаги 1),2) – прямой ход метода прогонки, 3),4) – обратный ход метода прогонки.

Теорема. Пусть коэффициенты a_i, b_i системы уравнений при i =2, 3, …, n–1 отличны от нуля и пусть

при i =1, 2, 3, …, n. Тогда прогонка корректна и устойчива.

При выполнении этих условий знаменатели в алгоритме метода прогонки не обращаются в нуль и, кроме того, погрешность вычислений, внесенная на каком либо шаге вычислений, не будет возрастать при переходе к следующим шагам. Данное условие есть ни что иное, как условие диагонального преобладания.

Для нашей краевой задачи имеем :

Тогда: , ,

Для нашей задачи условие устойчивости имеет вид:

.

Пусть

. (6)

Тогда

Пример. Найти решение задачи:

Выпишем разностную схему

Условие устойчивости примет вид

Возьмем .

Тогда

Или

Формулы прогонки были получены для СЛАУ (1):

Здесь x замены на u.

Следовательно,

Решим СЛАУ методом прогонки. Вычисления оформим в виде таблицы:

I	ai	ci	bi	fi	alfai	betai	ui
1		51	35	0.2	0.6863	-0.0039	0.4701
2	15	51	35	0.4	0.8598	-0.0113	0.6906
3	15	51	35	0.6	0.9186	-0.0202	0.8164
4	15	51	35	0.8	0.9403	-0.0296	0.9107
5	0	-1		1			1.0000

Порядок вычислений по формулам (4):

Ответ в столбце u_i.

Если забыли формулы, то их можно легко вывести. Главное запомнить основную формулу:

Прямой ход

Обратный ход

На практике часто граничные условия могут иметь более общий вид.

Рассмотрим следующую краевую задачу:

Найти решение ОДУ 2-го порядка

,

удовлетворяющую краевым условиям:

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

Аппроксимация:

В результате получим разностную схему:

Или

Контрольное задание. Лабораторная работа 5.

Численное решение дифференциальных уравнений

Задание.

1. Решить дифференциальное уравнение аналитически и численно указанными методами для двух значений шага интегрирования h=0.01; 0.001. Результаты расчета вывести на экран и распечатать в виде таблицы.

2. Построить графики функций y(x) (5 графиков).

Варианты уравнений и методов их решения приведены в таблице

Оформление результатов расчета

Таблица

х	Решения уравнения, у(x)
	Аналит	Численное
		метод 1		Метод 2
		h=0.01	h=0.001	h=0.01	h=0.001

Варианты уравнений и методов их решения

Таблица

Вар.	Вид уравнения	Метод	Вар.	Вид уравнения	Метод
1	у'=(xy2+x)/(y-x2y)	1,4	14	у'=cos(t)-y	3,5
2	у'=(1-2x)/y2	2,4	15	y'=exp(bx)-ay	1,4
3	у'=(1-x2)/xy	3,4	16	У'=-2y/(y2-6x)	2,4
4	у'=(y2-y)/x	1,5	17	у'=1/(2x-y2)	3,4
5	y'=(1+y)/(tg(x)	2,5	18	у'=sec(x)- y tg(x)	1,5
6	у'=exp(x)-1	3,5	19	y'=(exp(x)-y)/x	2,5
7	y'=y ln(y)/sin(x)	1,4	20	у'=1+y/(x(x+1))	3,5
8	у'=(1+y2)/(1+x2)	2,4	21	у'=(y+yx2-x2)/(x(1+x2))	1,4
9	у'=4x-2y	3,4	22	у'=cos(x-y)	2,4
10	у'=x exp(-x2)-2xy	1,5	23	у'=3x-2y+5	3,4
11	у'=2x-y	2,5	24	у'=sin(x)-y	1,5
12	у'=exp(-x)-2y	3,5	25	у'=exp(x)-y	2,5
13	у'=exp(-x)-2x	1,4	26	у'=exp(2x)-1	3,5

Примечание. Значение параметров a, b и начальные условия y|_x=x0=y₀ выбрать cамостоятельно.

Содержание отчета:

1. Название, цель работы и задание.

2. Математическое описание, алгоритм (структограмма) и текст программы.

3. Результаты расчета, пять графиков зависимости y(x) и выводы по работе.