Примеры численного решения ОДУ

1. Получить численное решение уравнения , при , . Построить график решения.

Решение. Переменной а присваиваем значение необходимого дифференциального уравнения

> a:=diff(y(x),x$2)+13*diff(y(x),x)+10=0;

Для получения численного решения также используется функция dsolve. Для того, чтобы указать на то, что требуется численное решение, необходимо дописать дополнительный параметр numeric. По умолчанию в системе Maple используется в этом случае метод Рунге-Кутта 4–5 порядка. Для указания на конкретный метод дописывается параметр следующего вида: method=<метод>. О методах, используемых в Maple можно узнать из справки системы.

Затем переменной присваиваем результат функции. Параметр range используется для определения интервала интегрирования уравнения.

> re:=dsolve({a,y(0)=10,D(y)(0)=4},y(x),numeric, range=0..1);

Результат, возвращаемый функцией dsolve с таким набором параметров – процедура с одним параметром, вычисляющая список значений (функции и производных) в точке, определяемой вносимым в параметр значением.

Для построения графика решения используется функция odeplot,которая позволяет построить график как значений функции, так и производных до порядка уравнения. После решения следующим параметром может указываться необязательный параметр – список переменных, по которым будет построен график. odeplot входит в состав пакета plots. Пакеты в Maple подключаются функцией with(<параметр>[,<функция>]), где параметром является пакет, а далее (не обязательно) можно писать список функций этого пакета, которые необходимо подключить к программе. Если список опустить, будут подключены все функции

> odeplot(re,[x,y(x)],color=blue);

Для получения графика производной в данном случае следует указать следующий список переменных: [x,D(y)(x)]. Можно также воспользоваться формой производной с использованием diff. Если в списке указать три параметра, то по ним будет построен трехмерный график. Если список будет опущен в описании, система сама выберет систему координат.

Лабораторная работа №1
Решение ОДУ различными методами

Задача: решить уравнение следующего вида: , где и при .

1. Решение методом неопределенных коэффициентов

Для начала отменим все определения и зададим функции p(x) и q(x), используя оператор функциональной зависимости.

> restart;

q:=x->cos(x);

p:=x->-tan(x);

Обозначим переменными ur[1] и ur[2] заданное уравнение и его однородную форму. Через xp и yp обозначим условие , то есть . Функция lhs(<уравнение>) возвращает левую часть уравнения. rhs – правую.

> ur[1]:=D(y)(x)+p(x)*y(x)=q(x);

ur[2]:=lhs(ur[1])=0;

xp:=0;

yp:=2;

Получим общее однородное решение уравнения, возьмем его правую часть. Заменим в ней постоянную _С1 произвольной переменной С, зависящей от x и присвоим полученное выражение переменной igr. Один из способов производить замену – это использование функции eval с двумя параметрами, где первый – выражение, а второй – список равенств, в соответствии с которыми будет произведена замена.

> rhs(dsolve(ur[2],y(x)));

igr:=eval(%,[_C1=C(x)]);

Затем возьмем производную полученного общего решения, для того, чтобы, заменив и в полном уравнении, определить С(х). Обозначим производную переменной igr1.

> igr1:=diff(%,x);

Произведем замену. Для нахождения С(х) решим дополнительное дифференциальное уравнение по С(х) с граничными условием .

> eval(ur[2],[y(x)=igr,D(y)(x)=igr1]);

dsolve({lhs(%)=rhs(ur[1]),C(xp)=yp},C(x));

Таким образом мы получили функцию С(х). Подставив ее в общее решение, мы получим полное решение уравнения:

> 'y(x)'=eval(igr,[C(x)=rhs(%)]);

2. Решение методом замены

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

> restart;

PDEtools[declare](y(x),u(x),v(x),prime=x):

q:=x->cos(x);

p:=x->-tan(x);

xp:=0;

yp:=2;

Обозначим переменной ur[1] заданное уравнение.

> ur[1]:=D(y)(x)+p(x)*y(x)=q(x);

Преобразуем уравнение к виду с производной. Произведем замену . С помощью функции collect (collect(<выражение>,<переменная>)) соберем коэффициенты при . Полученное уравнение обозначим переменной urav:

> convert(ur[1],diff);

eval(%,[y(x)=u(x)*v(x)]);

urav:=collect(%,v(x));

Функции и мы можем выбирать произвольным образом. Легко заметить, что уравнение заметно упрощается, если . Это условие обозначим как ur[1], а оставшееся выражение – ur[2].