Mathcad: решение краевой задачи для дифференциального уравнения 2 порядка
Краевая задача — это такое дифференциальное уравнение (или даже система дифференциальных уравнений) с заданными линейными соотношениями между значениями искомых функций на начале и конце заданного интервала.
Частые применения дифуров (обыкновенных дифференциальных уравнений, ОДУ) - решение задач вариационного исчисления, оптимального управления, механики жидкости и газа, баллистики, теории упругости и т.д.
Например, очень распространена на практике (и достаточно универсальна) такая постановка: для дифференциального уравнения второго порядка, имеющего вид
u'' + p(t)*u' + g(t)*u = f(t), t принадлежит интервалу [a,b]
поставлена краевая задача
k1*u(a) + k2*u'(a) = A
l1*u(b) + l2*u'(b) = B
Здесь u(t) - искомое решение, A, B - заданные краевые условия, определяющие значения искомой функции на концах интервала, p(t), g(t), f(t) - заданные функции коэффициентов, а штрихами, как обычно, обозначаем производные.
Если не вдаваться глубоко в теорию (она есть, например, тут, правда, обозначения другие), то такие задачи решаются методом конечных разностей, суть которого вот в чём:
1. Область непрерывного изменения аргумента t (интервал [a,b]) заменяется дискретным множеством точек (узлов). Количество этих узлов мы будем обозначать m, так что у нас выйдет дискретный аргумент
ti = a + τ*i, i=0,1,...,m, τ=(b-a)/m.
2. Искомая функция u(t) непрерывного аргумента t приближённо заменяется функцией дискретного аргумента ti на заданной сетке (а производные можно заменить конечными разностями, просто
по определению производной). Такая функция называется сеточной.
3. Исходное дифференциальное уравнение заменяется уравнением относительно сеточной функции. Это и есть разностная аппроксимация.
Решение уравнения в нашей постановке сводится к решению СЛАУ с трёхдиагональной матрицей, размерность которой будет (m+1)*(m+1), а все ненулевые элементы сосредоточатся на главной диагонали и двух прилегающих к ней диагоналях.
Было бы нерационально решать такую систему "в лоб", для неё есть экономичный метод прогонки. Ну а мы реализуем его и конечно-разностный метод в целом при помощи MathCAD.
В тестовых целях возьмём на просторах инета вот такую задачу, для которой известно решение:
постановка задачи
Введём в документ MathCAD все исходные данные задачи в виде переменных:
определение исходных данных
Это позволит вам, задав свои данные, решить другую задачу.
Реализуем метод прогонки решения краевой задачи для ОДУ 2 порядка, написав подпрограмму-функцию Progonka. Для простоты она будет пользоваться описанными выше переменными, а не получать кучу параметров.
реализация метода прогонки для решения краевой задачи (ОДУ 2 порядка)
Как и при решении задачи Коши, возвращаем матрицу из 2 столбцов, первый из которых - значения ti в узлах сетки, второй - вычисленные значения yi дискретной функции.
В учебных целях определим функции u(t) - точное решение задачи и Solution - формирование таблицы значений функции точного решения в тех же точках, в которых мы искали решение численным методом.
Останется вызвать обе функции и сравнить, что получилось, взяв норму разности векторов точного и приближённого решения. Оба этих действия на скрине:
сравнение решения конечно-разностным методом и точного решения
По малому значению погрешности видно, что прогонка работает хорошо, при увеличении значения m вдвое погрешность уменьшается для тестовой задачи более, чем вдвое, вывод о порядке аппроксимации метода
сделайте сами :)
В прилагаемом документе Mathcad есть также график, сюда не вывожу.
Решение краевой задачи стандартной функцией sbval можно попробовать тут
15.11.2014, 16:16 [29917 просмотров]