18. Дифференциальные уравнения в частных производных.

Дифференциальные уравнения в частных производных (ДУЧП, англоязычное сокращение Maple - pde - Partional Differential equation) в программе решаются сложнее, чем ОДУ (как и на бумаге). Уравнение задаётся оператором присвоения. При записи уравнения необходимо указать зависимость искомой функции от аргументов, например: f(x,y,t).

Ниже рассмотрим примеры двух алгоритмов решения ДУЧП. 1-й способ (которым часто пользуются и на бумаге) - редукция ДУЧП к системе ОДУ, что позволяет далее использовать алгоритмы, рассмотренные в п. 17. 2-й способ - при помощи специальных команд, предназначенных именно для решения ДУЧП. Более подробные справки находятся в Help по контекстным вызовам pde и pdsolve, или по общей схеме справки: Glossary/ Mathematics/Differential equation - и далее. Стандартные обозначения и операторы, введённые в предыдущих пунктах, ниже не объясняются.

>

18.1. Разделение переменных.

Рассмотрим в качестве примера решение 1-мерного волнового уравнения в однородной среде. Далее f(x,t) - искомая функция, которая может иметь различный физический смысл, v - скорость волны.

> pde1:=diff(f(x,t), x, x)-1/v^2*diff(f(x,t), t, t)=0;

Будем искать решение pde1 в виде произведения.

> subs(f(x,t)=X(x)*T(t), pde1): simplify(%): pde2:=expand((%)/X(x)/T(t));

Разнеся члены равенства в разные части, и учтя независимость переменных х и t, можем порознь приравнять обе части равенства к одной и той же константе - k^2 (далее введём обозначение w = k*v). Получим 2 ОДУ.

> ode1:=diff(X(x),`$`(x,2))=-k^2*X(x); diff(T(t),`$`(t,2))=-k^2*v^2*T(t): ode2:=subs(k^2=w^2/v^2,%);

Для их решения используем алгоритмы п. 17, выбрав НУ.

> dsolve( {ode1, X(0)=A1, D(X)(0)=0}, X(x)); X1:=A1*cos(k*x):

> dsolve({ode2,T(0)=0,D(T)(0)=B2*w},T(t)); T1:=B2*sin(w*t):

Итак, решение исходного ДУЧП в виде произведения есть:

> fxt:=X1*T1: combine(%); f:=subs(A1=2*C/B2, %);

Оператор combine представил решение суперпозицией волн, распространяющихся в противоположных направлениях. Константы A1 и B2 входят в решение одной комбинацией, и могут быть выражены одной новой константой С. Смысл константы k выясняется из условия периодичности решения с длиной волны L.

> X2:=A1*cos(k*(x+L)): eq1:=X2/A1=X1/A1; eq2:=k*(x+L)=k*x+2*Pi; K=solve(eq2, k);

Вообще имеем набор K = 2*Pi*n/L (n - целое) и соответствующих частот (см. параметр AllSolution, п.п. 9 и 10).

>

18.2. Решение командой pdsolve.

Обратимся к тому же примеру 1-мерного волнового уравнения, что и в п. 18.1 (см. pde1). Формальное применение команды решения даёт суперпозицию двух функций, отличающихся знаком при х, что обычно истолковывают как 2 волны, распространяющиеся в противоположных направлениях.

> pdsolve(pde1);

Задача состоит в отыскании этих двух функций конкретно. Для этого используем развёрнутую команду pdsolve с дополнительными параметрами (см. Help). Параметры INTEGRATE и build уточняют задание для программы. Параметр HINT (подсказка) предлагает программе искать решение в виде произведения (`*`), как и в п. 18.1.

> pdsolve(pde1, f(x,t), HINT= `*`, INTEGRATE, build); (1)

Выражение (1) содержит 4 константы интегрирования, заданный параметр v и неопределённый параметр, обозначенный программой _с1 (вообще, комплексный). Обозначим (1) (ниже ряд промежуточных результатов не выводятся на экран):

>F:=exp(_c[1]^(1/2)*x)*_C1*_C3*exp(_c[1]^(1/2)*v*t)+exp(_c[1]^(1/2)*x)*_C1*_C4/exp(_c[1]^(1/2)*v*t)+1/exp(_c[1]^(1/2)*x)*_C2*_C3*exp(_c[1]^(1/2)*v*t)+ 1/exp(_c[1]^(1/2)*x)*_C2*_C4/exp(_c[1]^(1/2)*v*t):

Константы _С1 -_С4 найдём из Н.У., которые наложим на F и её производные по x и t. Ниже dx и dt - не дифференциалы, а просто обозначения!

> Fdx:=diff(F,x): Fdt:=diff(F,t):

Эти условия зададим несколько по-иному, чем в п. 18.1, что упрощает запись, но не имеет принципиального значения. Именно, потребуем, чтобы при (t=0, x=0) искомая функция равнялась А, а обе производных - нулю:

> subs([x=0, t=0],F): subs(exp(0)=1, %): eq1:=%=A;

> subs([x=0, t=0],Fdx): subs(exp(0)=1, %): simplify(%): eq2:=%=0;

> subs([x=0, t=0],Fdt): subs(exp(0)=1, %): simplify(%): eq3:=%=0;

(Для упрощения сделана подстановка exp(0)=1, что программа сама не делает, но знает). Получили систему 3-х уравнений, позволяющих выразить три константы через одну, напр. _С1 (см. п. 9), что и подставим в F. (Надо точно сохранять введённые программой обозначения констант!)

> solve({eq1, eq2, eq3}, {_C2,_C3,_C4});

>subs([_C2=_C1,_C3=1/4*A/_C1,_C4=1/4*A/_C1],F):simplify(%):expand(%):combine(%); (2)

В итоге _C1 выпало из решения, вместо него вошло заданное А. Решение (2) описывает весьма широкий класс функций, отличающихся параметром _с1. Будем искать периодические решения, для чего положим _с1<0, так что корень из этой величины - мнимый.

> factor(%);

> subs(_c[1]^(1/2)=I*k,%); convert(%, trig);

> y:=simplify(1/4*A*(2*cos(k*(v*t+x))+2*cos(k*(v*t-x))));

Полученное решение согласно с решением, данным программой в ответе на самую первую команду, и отличается от решения, найденного в п. 18.1 только сдвигом по фазе. Дополнительное определение k - как в п. 18.1.

Решение ДУЧП с заданными граничными условиями здесь не обсуждается.

>