MathCAD: решаем дифференциальное уравнение 1 порядка

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

Итак, дано уравнение вида dy/dx=D(x,y), интервал поиска решения x∈[a,b] с начальным условием y(а)=const. Интервал будем разбивать с шагом h=(b-a)/n, где n - число значений, которые принимает дискретизованная переменная x.

Сначала определим функцию, задающую производную, начальное условие, интервал поиска решения и количество узлов сетки - эти данные понадобятся всем методам. Так как задача - учебная, добавим и определение функции точного решения:

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

Метод Эйлера легче всего запрограммировать, ведь он состоит в простейшем линейном приближении производной: из исходного уравнения, по определению производной, получаем (y[i+1]-y[i])/(x[i+1]-x[i])=D(x,y), откуда, при постоянном шаге по x, равном h=x[i+1]-x[i], имеем y[i+1]=y[i]+h*D(x,y).

Метод Эйлера в MathCAD

Метод возвращает то же, что стандартная функция для решения дифференциального уравнения (см. ниже) - матрицу из 2 столбцов, первый из которых - значения x[i] в узлах сетки, а второй - вычисленные значения y[i].

Для проверки полученного решения оценим норму разности между вычисленными подпрограммой значениями y[i] и значениями, полученными функцией точного решения f(x) в тех же значениях аргумента x[i]. Также построим график.

Иллюстрация решения методом Эйлера

Этот метод имеет первый порядок точности, о методике практической оценки порядков точности см. здесь.

Модифицированный метод Эйлера, он же метод "предиктор-корректор", он же метод Эйлера-Коши, имеет второй порядок точности и умеет уточнять решение, полученное обычным методом Эйлера, за счёт дополнительного вычисления значений D(x,y) "между" узлами сетки:

Модифицированный метод Эйлера в MathCAD

Наконец, наиболее популярный на практике метод Рунге-Кутты 4 порядка реализован в стандартной функции MathCAD, воспользуемся ей вместо "ручного" программирования:

Метод Рунге-Кутты 4 порядка в MathCAD

Не сложнее было бы запрограммировать метод Рунге-Кутты 4 порядка и "вручную":

Метод Рунге-Кутты 4 порядка в MathCAD (программа)

В приложенном документе видно, что все результаты совпали.

Увеличивая n, например, до 20 и 40, можно проверить, соблюдаются ли теоретические порядки точности методов.

 Скачать документ MathCAD 15 (.xmcd) в архиве .zip с этими расчётами (50 Кб)

P.S. По просьбам трудящихся: как проверить, является ли известная функция u(x) решением дифференциального уравнения dy/dx=D(x,y) ?

В нашем случае, так как функции D(x,y) и u(x) уже определены в документе, зададим и символьно оценим функцию первой производной от u(x), обозначим её, например, z(x).

Затем перепишем уравнение в виде z(x)-D(x,u(x))=0 (подставили вместо y возможное решение u(x) и перенесли D(x,y) в левую часть уравнения) и символьно оценим его оператором simplify с панели инструментов "Символика".

Если результат оценивания равен нулю, уравнение превратилось в тождество и u(x) является его решением, для "неправильной" же u(x) нуля в операторе simplify не получится.

Можно было и сразу записать производную в операторе simplify, воспользовавшись значком производной с панели инструментов "Математический анализ". Оба способа показаны на этой картинке:

Mathcad, как проверить, является ли функция решением дифференциального уравнения

17.11.2013, 21:57 [62916 просмотров]