Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 100
Текст из файла (страница 100)
Решение ОДУ третьего порядка (рис. 14.4)Given'"(t) + y"(t) + y'(t) + y(t) = -t-sin(t)y"(0) = 2y'(0)=ly(0)=0у :=odesolve(t,30,30000)200 т~ . /\10y(t)-200 -1tРис. 14.4. Модель резонирующей колебательной системыИзучая пример 14.8, вы, наверное, заметили, что дифференциальное уравнение в немзадано иначе, чем в примере 14.7. Вообще же, обе эти формы абсолютно эквивалентны,14.2.
Численное решение ОДУ в форме задачи Коши * 4 4 7и выбор одной из них должен определяться прежде всего вашими предпочтениямии той областью, к которой относится решаемое уравнение.Важным условием успеха при численном решении дифференциальных уравнений является правильный выбор величины шага. Так, если он будет определен недостаточномалым, то найденное решение может быть весьма и весьма далеким от истинного. Приведем пример ошибочного решения ОДУ при слишком большой ширине шага.Пример 14.9.
Осцилляция при недостаточно малой длине шага (рис. 14.5)Giveny'(t) = l -у(5) = 5y:=odesolve(t,50,10)yl:=odesolve(t,50,1000)'82-Ш т.0000000520tРис. 14.5. Ошибочное (слева) и верное (справа) решение дифференциального уравненияКонечно, в приведенном примере возникшая ошибка (см. рис. 14.5) связана с искусственно завышенной длиной шага, и даже при использовании его величины, принятойпо умолчанию (step=0.1), решение было бы найдено верно.На практике бывает необходимо решать такие дифференциальные уравнения, нужноезначение шага для которых может быть совсем не очевидно (и step=0.1 для них будетпорождать еще более значительные колебания или отклонения, чем step=4.5 в примере 14.9).
В тех случаях, когда у вас возникают сомнения в верности заданной вами длине шага, вы можете руководствоваться следующим простым правилом: ее можно принять как корректную, если уменьшение ее величины в 10 раз не приводит к изменениюрезультата в некоторой точке в пределах количества десятичных знаков, отвечающегонужному уровню точности.Как уже было отмечено выше, по умолчанию вычислительный блок Given-Odesolve использует популярный метод Рунге-Кутта 4-го порядка с постоянным шагом. ОднакоMathcad позволяет сменить его при необходимости на так называемый адаптивныйалгоритм (adaptive) или алгоритм решения жестких систем (stiff).
Чтобы это сделать,щелкните правой кнопкой мыши на функции odesolve. В открывшемся при этом контекстном меню переместите флажок из строки Fixed (Фиксированный) в строку Adaptive(Адаптивный) или Stiff (Жесткий). Отличие адаптивного метода от простого алгоритма Рунге-Кутта состоит в том, что используемая им длина шага не постоянна, а зависитот скорости изменения функции результата (что определяется, в частности, величиной4 4 8 • Глава 14.
Дифференциальные уравненияпервой производной). Такой подход может быть более эффективен в случае не слишком хорошо подобранного пользователем шага, а также для решения сложных дифференциальных уравнений с быстро изменяющимися и осциллирующими функциями решений. Если же вам нужно решить жесткую систему уравнений или систему,содержащую алгебраические ограничения, следует воспользоваться методом Stiff(рис. 14.6).у := odesolv)*'* W 1 rh IFixedРис. 14.6. Смена используемого численного метода14.2.2. Решение ОДУ с помощью встроенных функцийВ более старых версиях Mathcad для решения ОДУ использовались специальныевстроенные функции. Однако с появлением более наглядного и простого способа, связанного с применением вычислительного блока Given-Odesolve, они уже не представляют особой практической важности.
Несмотря на это, о возможности их применениядля решения одинарного ОДУ нужно иметь четкое представление, так как в некоторых случаях (например, в программировании) использовать вычислительный блокнельзя.Так как для решения одного дифференциального уравнения соответствующие встроенные функции применяются редко, то подробно описывать их в этом разделе мы небудем и сразу приведем соответствующий пример (в случае необходимости читательбез труда в нем разберется). Обстоятельно же об указанных функциях мы поговоримв разделе, посвященном решению систем ОДУ.Пример 14.10.
Альтернативный метод решения ОДУ (рис. 14.7)уО:=1М:=100D(t,y):=sin(t)-yxl:=0yl :=rkfixed(yO,xl,x2,M,D)х2:=10y2:=Bulstoer(yO,xl,x2,M,D)уЗ :=Rkadapt(yO,xl,x2,M,D)14.2. Численное решение ОДУ в форме задачи Коши * 4 4 910yl=012300.10.20.310101.005 у2 = 11.0221.046300.10.20.31011.005 уЗ=1.021.046012300.10.20.311.0051.021.046Рис. 14.7. График решения ОДУ, полученного с помощью встроенной функции rkfixedПрименять приведенные встроенные функции можно лишь для решения дифференциальных уравнений первого порядка, как линейных, так и нелинейных. Найти же с ихпомощью решение уравнения более высокого порядка, в отличие от использования блока Given-Odesolve, непосредственным образом невозможно. Впрочем, это можно сделать,сведя указанное уравнение к системе дифференциальных уравнений. Приведем пример решения ОДУ второго порядка с помощью вычислительного блока и встроеннойфункции rkfixed (реализующей метод Рунге-Кутта четвертого порядка с фиксированным шагом).Пример 1 4 .
1 1 . Альтернативные способы решения ОДУ второго порядка(рис. 14.8)Вычислительный блок:GivenУ(0)=1У'(О)=3z:=odesolve(t,5,100)Сведение к системе ОДУ:у:=D(t,y):=у + 2уr:=rkfixed(y,0,5,100,D)04 5 0 • Глава 14. Дифференциальные уравненияtг<о>Рис. 14.8. Альтернативные решения ОДУ, полученные с помощью вычислительного блокаи встроенной функцииИзучив пример 14.11, вы, наверное, согласитесь, что решать дифференциальные уравнения, особенно высокого порядка, гораздо проще и лучше с помощью вычислительного блока, чем встроенных функций. Действительно, зачем выполнять лишнюю, иногда весьма сложную и объемную, работу по разложению уравнения в систему ОДУ, еслиэто отлично может сделать и сама программа.
К тому же решение при использованиивычислительного блока получается гораздо более наглядным и доступным пониманиюдля не владеющего Mathcad человека, чем при применении системных функций с ихзапутанным синтаксисом. Поэтому, на мой взгляд, прибегать к встроенным функциямдля решения ОДУ стоит лишь в том случае, если использовать блок Given-Odesolve невозможно.14.2.3. Системы ОДУНачиная с 11 версии Mathcad возможности вычислительного блока Given-Odesolve значительно расширились.
Теперь с его помощью можно находить численные решения нетолько ОДУ, но и их систем. Как и к ОДУ, к уравнениям системы предъявляются следующие требования: каждое из них должно быть линейным (высшая производная недолжна содержать каких-либо сомножителей) и иметь соответствующее количествоначальных или граничных условий.К сожалению, вычислительный блок не может решать системы нелинейных дифференциальных уравнений.
Однако в Mathcad имеются специальные встроенные функции,позволяющие находить решения как линейных, так и нелинейных систем ОДУ. Всеготаких функций три.• Rkfixed(yO,tO,tl,M,D). Задействовав эту системную функцию, можно решить задачуКоши с помощью, пожалуй, самого популярного из численных алгоритмов: методаРунге-Кутта 4-го порядка с фиксированным шагом.
Эта функция подходит длякачественного и быстрого решения подавляющего большинства систем ОДУ.• Rkadapt(yO,tO,tl,M,D). Эта функция реализует адаптивный алгоритм Рунге-Кутта. Еепринципиальным отличием от rkfixed является то, что при вычислении соответствующих приближений она использует не постоянный, а зависящий от скорости изменения функций решения шаг дискретизации переменной. Однако, несмотря на это,в качестве ответа Rkadapt возвращает значения функций в равномерно распределенных, исходя из заданной пользователем величины промежутка, точках интервала.Использовать же рассматриваемую функцию следует в случае более жестких и слож-14.2.
Численное решение ОДУ в форме задачи Коши*451ных систем уравнений — это позволит повысить точность расчета или сэкономитьвремя по сравнению с применением простого метода Рунге-Кутта.• Bulstoer(yO,tO,tl,M,D). Метод Булирша-Штера. Использовать этот алгоритм стоитв том случае, если вы уверены, что функции решения вашей системы достаточногладкие и плавно изменяющиеся. При выполнении этого условия функция Bulstoerпозволяет получать более точные решения, чем rkfixed, затрачивая на это меньшевремени.Для корректного использования описанных выше функций система дифференциальных уравнений должна быть записана в векторном виде:= D(Y(x),x)где Y-(x) — вектор первых производных системы, D(Y(x),x) — вектор-функция, каждая строка которой содержит левую часть соответствующего уравнения системы.Кроме того, в векторной форме должны быть определены начальные условия:Разобравшись с принципом приведения системы ОДУ к векторной форме, вы без труда сможете задать и параметры для рассматриваемых встроенных функций:•уО — вектор начальных условий; в нем вы должны определить числовые значенияискомых функций на левой границе'интервала изменения переменной;• tO — начальная точка для переменной;• tl — конечная точка расчета;• М — количество шагов, при котором численный метод будет решать систему;D — вектор-функция, содержащая левые части уравнений системы.
Должна бытьзадана как функция двух переменных: скаляра t и вектора у (то есть все искомыефункции системы должны быть представлены как элементы одного вектора).Результатом работы приведенных функций является матрица, в первом столбце которой содержатся узловые величины переменной t, а в остальных — значения неизвестных функций системы, рассчитанные в этих точках. При этом порядок расположениястолбцов с найденными величинами искомых функций определяется последовательностью, в которой они были занесены в вектор у.•Аппарат дифференциальных уравнений широко используется для описания не толькофизических, химических или биологических процессов, но и различных явлений в медицине, экономике, демографии и множестве других современных наук. Зачастую припопытке решения той или иной системы уравнений можно получить совершенно неожиданный результат. В этом разделе мы рассмотрим наиболее интересные примерырешения систем ОДУ: экологическую модель «хищник — жертва», некоторые видыдинамических систем, а также движение ракеты в поле тяготения небесных тел.Пример 14.12.
Модель «хищник — жертва» (Лотки-Вольтерра)Модель «хищник — жертва» была предложена независимо друг от друга американским физикомАльфредом Лоткой в 1925 году и итальянским математиком Вито Вольтерра в 1926 году. Онаописывает эволюцию численности взаимодействующих популяций хищников и жертв на протяжении определенного промежутка времени. Рассмотрим количественные изменения, происходящие4 5 2 • Глава 14. Дифференциальные уравненияв популяциях рысей и зайцев. Зайцы, питаясь растительностью, размножаются с постоянной скоростью А, в результате чего их численность возрастает:— x(t)=AxdtРыси поедают зайцев, что уменьшает их численность, причем скорость этого процесса By пропорциональна количеству хищников у:—x(t)=-BxydtОбщая динамика популяции зайцев описывается уравнением:—x(t) = A x - ВхуdtВ результате количество зайцев становится настолько большим, что приводит к резкому увеличению количества рысей.
Популяция хищника растет пропорционально имеющейся добыче:- y ( t ) = BxydtУмирают рыси естественной смертью или под влиянием внешних факторов:-y(t)=-cydtВ итоге количество хищников определяется уравнением—y(t) = - C y + ВхуdtНа определенном этапе рысей становится так много, что количество зайцев быстро уменьшаетсявследствие интенсивного поедания хищниками. Резкое сокращение запасов пищи вызывает снижение численности рысей, что возобновляет рост популяции зайцев.
В конечном счете, процессповторяется заново. Ниже приведено решение и фазовые портреты полученной системы уравнений для двух начальных условий (рис. 14.9).B:=lA:=1.5А-G(t,y):=Byyo- - o- iу0:=Вуо-у,-С)С:=310yl:=r2:=rkfixed(yl,0,100,5000, G)rl :=rkfixed(yO, 0,100,5000,G)1020010510.020.040.060.080.19.2588.457.6066.7575.9355.716.4197.0997.7188.251234514.2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
