Часто Maple не сможет дать решение ДУ, выраженное в известных вам функциях. В этом случае потребуется численное решение и, возможно, его график либо решение при каком-то значении переменной, например при x = 2, и т. п. Для этого применяется команда dsolve с опцией type=numeric, как показано в примере ниже. Но прежде чем начать действия, следует предупредить, что:

(1) «Численное» подразумевает, что применяются числа с плавающей запятой. Если для начальных условий задать целые числа или определения констант, то dsolve(....type=numeric) может не сработать.

(2) Требование получения чисел не учитывает «нелюбовь» Maple к такому способу вычислений. Поэтому возможны трудности с выполнением полезных действий над результатом.

(3) Будет намного лучше, если задавать Maple ДУ 1-го порядка, а не 2-го, 3-го.

(4) Результат выполнения dsolve присваивайте переменной, чтобы можно было сослаться на него в дальнейшем: F:=dsolve(...type=numeric).

Начнем с численного решения для гармонического осциллятора. Вначале определим уравнения для x(t) и v(t):

• restart:omega:=1.;
• Eqx:=diff(x(t),t)=v(t);Eqv:=diff(v(t),t)=-omega^2*x(t);

затем определим начальные условия:

• IC:=x(0)=1.,v(0)=0.;

и в конце вызовем dsolve с численной опцией:

• XV:=dsolve({Eqx,Eqv,IC},{x(t),v(t)},type=numeric);

Теперь XV – это процедура, которая ожидает, чтобы ей дали значение t, а она вернет приближенные решения для x(t) и v(t):

• XV(2.);

Учтите, что ответ дается с точностью 18 десятичных знаков. В нашем случае он известен заранее (x(t) = cos(t)), поэтому возможна проверка (после увеличения числа знаков Maple останется в режиме с плавающей запятой).

• Digits:=18:cos(2.);

Очевидно, Maple врет – результат dsolve вполне точен (7 значащих цифр), но не настолько точен, как должны бы показать возвращаемые им числа.

Увеличим точность, указав, какой численный метод следует применить, запросив справку ?dsolve, numeric. Хороший выбор – rkf45. Если не подходит, выберите другой и задайте допустимую погрешность, как ниже:

• XV:=dsolve({Eqx,Eqv,IC},{x(t),v(t)},type=numeric,method=rkf45,abserr=1e-10);
• XV(2);

Это лучше, но Maple не всегда это позволяет. Попробуйте в команде выше изменить величину погрешности на 10-15 и посмотрите, что получится. (Не оставляйте ее столь малой, загоните комп.)

Совет: не пытайтесь вычислять со слишком высокой точностью, например 1 часть на 1 000 000 000 000. Обычно хороший результат получается при 5–6 значащих цифрах. Чтобы удостовериться, что примененные параметры хороши для решаемой задачи, выполните численный расчет дважды с разной погрешностью.

Созданная процедура XV решает ДУ численно. Теперь построим график решения. Здесь нюанс: assign(XV) не сработает, как было в символьных вычислениях, т. е. эту конструкцию нельзя использовать. Проще всего применить команду odeplot, которой надо передать результат, полученный из команды dsolve(....type=numeric). Но вначале напишите: with(plots).

Прежде чем продолжить, обсудим еще одну проблему.

В команде dsolve установлена погрешность 10-10, что слишком мало и при запросе к odeplot потребует много точек. Результатов придется долго ждать. Поэтому надо вначале установить приемлемую погрешность для рисования графика:

• XV:=dsolve({Eqx,Eqv,IC},{x(t),v(t)},type=numeric,method=rkf45,abserr=1e-6);

Чтобы нарисовать x(t), применим:

• with(plots):odeplot(XV,[t,x(t)],0..20,numpoints=500,color=navy);

Для рисования x(t) и v(t):

• dsolve odeplot(XV,[[t,x(t)],[t,v(t)]],0..20,numpoints=200);

или создать фазовый график, т. е. v(t) в зависимости от x(t) со временем t в качестве параметра:

• odeplot(XV,[x(t),v(t)],0..20,numpoints=200);

Отметьте себе:

(1) в odeplot есть все обычные опции plot, так что можно улучшить вид графика;

(2) если ошибка задается малой, а numpoints задано большим, то придется долго ждать построения графика.

Раньше эту задачу не удалось решить символьно. Теперь сделаем численный расчет и построим y(x) от x = 0 до x = 20.

В Maple есть пакет DEtools со множеством возможностей для решения ДУ, одна из которых визуально показывает путь решения ДУ («потоки»). Он называется DEplot. Ниже показан пример применения (опять осциллятор с затуханием):

где – угловая частота, g – скорость затухания.

Команды ниже показывают, как включить набор связанных ДУ в команду DEplot, чтобы показать решение в виде кривой в [x,v] плоскости, называемой в физике фазовым пространством. С ним часто имеют дело в классической механике.

Вызовем пакет DEtools:

• restart:with(DEtools):
• eqx:=diff(x(t),t)=v(t);eqv:=diff(v(t),t)=-x(t)-.5*v(t);

Теперь решим ДУ численно с начальными условиями x(0) = 1 и v(0) = 0. Запустим на 20 секунд и построим решение в фазовом пространстве [x,v]. DEplot будет рисовать стрелки, которые показывают ход решения. Это красивое свойство позволяет визуализировать решение при заданных начальных условиях.

DEplot([eqx,eqv],[x(t),v(t)],t=0..20,[[x(0)=1,v(0)=0]],
x=-1..1,v=-1..1,stepsize=.01,scaling=constrained,color=red,
linecolor=navy);

Наберите запрос ?DEplot, чтобы изучить команду и понять, что делает каждый ее член.

Затем увеличьте скорость затухания и проследите за изменением картинки. И в конце замените член затухания с gv на gv|v| (затухание в воздухе) и снова посмотрите на картинку. Объясните изменения.

Бывает, что требуется искать не решение ДУ с известными начальными условиями, а нужен его неизвестный параметр, значение которого определяется из решения в другой точке. Или известны только некоторые начальные условия, и надо подгонять неопределенные начальные условия до тех пор, пока решение не удовлетворит конечным условиям. Подгонка такого рода напоминает артиллерийскую пристрелку: стреляем из точки а и попадаем в мишень в точке b (или рядом с ней). Если получается символьное решение, то подогнать параметры несложно, но для численного решения это весьма непросто, и возможностей Maple может не хватить. Ниже расскажем о процедуре, показывающей пути решения проблемы.

Вначале сформулируем конкретную задачу.

Уравнение Шредингера (УШ) имеет вид:

Переведем в более приятный безразмерный вид, заменив x на ξ:

x = aξ,  где  ,

и, заменив E на собственное значение λ:

УШ в новых переменных:

Решения этого уравнения бывают двух типов: четные либо нечетные функции от ξ. Следовательно, при ξ = 0 есть два разных начальных условия:

четное –        ψ(0) = 1 и

нечетное –    ψ(0) = 0 и

Конкретные значения параметров не очень важны. Важнее суть решения. Поэтому нет ничего страшного, если, к примеру, начальный наклон кривой задан =1, а не 10-16. В ДУ видно, что в каждое слагаемое входит ψ, поэтому если решение найдено, то удвоенное значение тоже будет решением. Поскольку «размер не важен», берем 1. В реальных задачах волновая функция нормируется, и тогда все хорошо. (Узнаете из квантовой механики.)

Решение для волновой функции (в.ф.) имеет физический смысл, если только в.ф. стремится к 0 при ξ→∞. Такое возможно только для совершенно определенных значений параметра λ (собственных значений). Поэтому нужна процедура, которая:

(1) умеет решать ДУ численно и

(2) умеет подгонять параметр λ так, чтобы в.ф. ψ была мала при больших ξ.

Вот Maple-овский код, содержащий процедуру пристрелки для решения более простой похожей задачи (ДУ для гармонического осциллятора):

с y(0) = 0 и

Нужно подогнать λ, применив условие y(1) = 0. Запустите процедуру (ниже), затем модифицируйте ее, чтобы в УШ найти значения λ для первых шести связанных состояний (3 четных и 3 нечетных). Придется решать проблему, насколько далеко надо продолжать («тянуть») ξ, чтобы аппроксимировать бесконечность. Не переоцените! Иногда достаточно 5–10.

Начнем заново и загрузим odeplot:

• restart:with(plots):

Определим ДУ и присвоим их deqn. (В этом разделе кода варьируемый параметр, удовлетворяющий условиям «стрельбы», назван lambda.)

• deqn:=diff(y(x),x$2)+lambda^2*y(x)=0;

Определим левую и правую конечные точки интервала:

• a:=0;b:=1;

Определим начальные условия в начале интервала:

• bc:=y(a)=0,D(y)(a)=1;
• yright:=0;
• lambda:=2.5;

Учтите, что опция startinit=true говорит dsolve всегда стартовать с определенных начальных условий. При использовании fsolve важно иметь набор, чтобы найти значение параметра, который даст решение при y(b) = yright.

• S:=dsolve({deqn,bc},y(x),type=numeric,startinit=true);
• odeplot(S,[x,y(x)],a..b);
• F:=proc(z)
• local a,b,deqn,bc,S;
• deqn:=diff(y(x),x$2)+z^2*y(x)=0;
• a:=0;b:=1;
• bc:=y(a)=0,D(y)(a)=1;
• if type(z,numeric) then
• S:=dsolve({deqn,bc},y(x),type=numeric,startinit=true,
  abserr=1e-6);
• subs(S(b),y(x))
• else
• 'F(z)';
• end if;
• end;

Теперь можно перенести F(λ) в fsolve, чтобы найти значение λ, которое удовлетворяет F(b)=yright. Для этого установите Digits в наименьшее значение, чтобы fsolve не перетрудился, пока ищет правильный ответ. При этом помните, что dsolve не всегда удовлетворяет требованию большей точности. Если вы не в процедуре, найдите правильный баланс между abserr в dsolve и поставленными здесь установками Digits, иначе fsolve будет работать бесконечно.

• Digits:=5;
• F(3.);
• lambda:=fsolve(F(z)=yright,z=3.0..3.2);
• S:=dsolve({deqn,bc},y(x),type=numeric):
• odeplot(S,[x,y(x)],a..b);

Теперь решение на правом конце стремится к нулю, как и хотелось.