Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 103
Текст из файла (страница 103)
Согласно закону Кеплера, она также будет двигаться по одному из конических сечений — параболе или гиперболе.Решение ищется аналогичным образом (рис. 14.16). Изменим лишь составляющие начальнойскорости.Vx:=lVy:=2Vz:=3(х,у ,z) ,CresteSpace(F ,0, lJO ,200фРис. 14.16. Траектория движения ракеты в поле тяготения одной планеты при ненулевойначальной скорости, направленной не к центру планетыПерейдем к более сложному случаю, когда ракета, имеющая нулевую начальную скорость, движется в поле гравитационных сил двух планет.Особенностью такого движения является то, что траектория ракеты всегда будет находитьсяв одной плоскости с центрами тяжести планет независимо от их массы и расположения.Изменим входные условия для масс и координат. Остальные выражения останутся неизменными.М-=|У:=.*•=.00001-юз07012-1031402080010120Решение системы (рис.
14.17) полностью подтверждает сказанное выше.И, наконец, рассмотрим последний случай, иллюстрирующий путь ракеты с ненулевой начальной скоростью в поле сил притяжения пяти планет.4 6 4 • Глава 14. Дифференциальные уравнения(x,y,z),CreateSpace(F,O,100,2000)(x,y,z),CreateSpace(F,O,100,2000)Рис. 14.17. Движение ракеты с нулевой начальной скоростью в поле сил гравитации двух планетВведем массы и координаты планет, а также начальную скорость ракеты. В маркере функцииOdesolve укажем конечную точку интегрирования 2500. Остальные параметры и уравненияоставим без изменения.М:=001234у.=01-1032-1033-юз2-1031-10301234Vx:=l900700500300100Vy:=2z:=02006001004005001234001234100200500300100Vz:=3Траектория ракеты будет представлять собой сложную кривую (рис.
14.18).(x,y,z) ,CreateSpace(F,0,2300,200G)O,y,z), CreateSp ac e(F, 0,2300,2000)Рис. 14.18. Путь ракеты в поле сил притяжения пяти планет (начальная скорость не равна нулю)Заканчивая разговор о решении систем ОДУ в Mathcad, хотелось бы посоветовать: относитесь к результатам работы встроенных функций и вычислительного блока критичнои осторожно. Всегда старайтесь искать решение с помощью двух различных алгоритмов и принимайте результат как истинный, только если их ответы совпадут.
Параметры встроенных функций настраивайте так, чтобы их изменение в 2-10 раз не приводи-14.2. Численное решение ОДУ в форме задачи Коши.;. 4 6 5ло к изменению в ответе. Руководствуясь описанным подходом, вы сможете избежатьмногих неприятных сюрпризов. В следующем разделе мы подробно разберем численные методы, лежащие в основе работы описанных выше встроенных функций.14.2.4.
Численные методы, применяемые встроеннымифункциями для решения ОДУ и их системМетоды Рунге-КуттаМетод Рунге-Кутта 4-го порядка наиболее популярен при решении задачи Коши, поскольку является достаточно точным и не требует вычисления производных высшихпорядков, в отличие, например, от метода Тейлора, что значительно сокращает времярасчетов. Чтобы понять, каким образом функция rkfixed, реализующая метод, находитчисленное решение, и насколько оно оказывается верным, попытаемся проанализировать сам алгоритм. Очевидно, самый лучший способ сделать это — самостоятельно написать соответствующую программу.= 2Рассмотрим уравнение y'(t) t ~y, для которого имеется аналитическое решение, чтобы сравнить, насколько близки к истинным значения, полученные численным методом.
Зададим начальные условия y(to)=yO и интервал поиска [t o ;tj, который разобьемна множество подынтервалов [t;t. +1 ]. Их количество определяется величиной выбранного шага h. Именно на этих промежутках мы будем вести поиск интегральной кривой,описывающей функцию решения. Сущность его заключается в нахождении семействатангенсов угла наклона кривой на данном промежутке, которое затем используется дляприближенного расчета приращения функции.
Сначала определим направление движения на заданном подынтервале. Для этого необходимо найти наклон интегральнойкривой в точке t . Очевидно, что он будет равен тангенсу угла наклона касательнойк нашей кривой в этой точке:Затем, двигаясь по касательной, сделаем половинный шаг и попадем в точку с абсциссой t.+h/2. Понятно, что приращение искомой функции составит h/2-kl.
Теперь определим наклон кривой в новой точкек2 = ff t . + - , y . + - - k lи вновь сделаем половинный шаг из исходной точки, правда, в направлении, заданном к2.Абсцисса останется той же, а ордината примет вид y.+h/2-k2. Повторим уже знакомуюнам операцию: найдем тангенс угла наклона касательной к функции в этой точке:k3 = ff t. + —, у. Н— • к2который используем для определения наклона функции на правой границе интерваларазбиения:h- k3)Итак, мы получили семейство тангенсов угла наклона кривой в трех точках интервала:kl — в начальной точке, к2 и кЗ — в середине и к4 — в конце.
Эти значения необходимы для того, чтобы провести численное интегрирование функции, задающей наклон4 6 6 •:• Глава 14. Дифференциальные уравненияна промежутке [t^t^J методом Симпсона, о котором мы говорили ранее. Функцияв этом случае интерполируется параболой. Чтобы определить ее вид, достаточно трехточек подынтервала, в которых значение функции известно. Выше мы определиличетыре таких точки: kl, k2, k3 и к4, поэтому в середине промежутка наклон кривой зададим как среднее арифметическое к2 и кЗ. Применив формулу Симпсона с половинными шагом, мы найдем приращение искомой функции на данном подынтервале.Эта идея лежит в основе метода Рунге-Кутта 4-го порядка, который позволяет подобным образом рассчитать приращение функции на всем интервале [t o ;tj.
Исходя из условия y(to)=yO, мы нашли решение yi в точке t . Аналогично, приняв за начальное условие y(t)=y j , можно определить функцию в точке t.+1 и т. д. В общем виде- • (kl 2k2 + 2k3 + k4)Чтобы наглядно представить всю последовательность действий, проведенных намивыше, напишем программу, реализующую данный метод.Пример 14.18. Программная реализация метода Рунге-Кутта 4-го порядка(рис.
14.19)f(t,y):=t2-yDif(f,a,b,yO,M):for ieO.. Mkl4-f(t,y.)h2h> У1+2t + - , y . + --k22 > 2(kl+2-k2+2-k3+ k4)^. «- t t <r- t + h)rуC:=Dif(f,0,10,1,100)iif14.2. Численное решение ОДУ в форме задачи Коши• 467Рис. 14.19. Решение ОДУ с помощью приведенного алгоритмаВ корректности работы программы легко убедиться, сравнив решения, полученныес помощью встроенных функций и нашего алгоритма.На каждом шаге алгоритм аппроксимирует дифференциальное уравнение по приведеннойсхеме, поэтому попытаемся проанализировать, какие факторы влияют на точность полученных с его помощью результатов.
Очевидно, возникает ошибка приближения искомой функции параболой при использовании формулы Симпсона, а также погрешностьусреднения значений к2 и кЗ в центральной точке подынтервала. Такой подход к оценкепогрешности можно назвать геометрическим. С аналитической точки зрения метод Рунге-Кутта 4-го порядка использует для нахождения решения только первые пять слагаемых, полученных при разложении искомой функции в ряд Тейлора в окрестностизаданной точки. Если на отдельном шаге погрешность пренебрежимо мала, то на интервале происходит накопление ошибки, что отражается на точности результата.
Поэтому важно выбрать оптимальное соотношение между длиной шага и их количеством, при котором погрешность аппроксимации окажется минимальной, а время расчета приемлемым.Как правило, с уменьшением длины шага точность возрастает до определенной величины, которая не меняется при дальнейшем разбиении интервала. Проиллюстрируем сказанное. В табл. 141 приведены истинные решения задачи Коши y'(t)=t2—у, у(0)=1 наинтервале [0;10], а также рассчитанные нашим алгоритмом при различной длине шага.Таблица 1 4 . 1 .
Решения задачи Коши y'(t)=t2—у, у(0)=1 на интервале [0; 10] методомРунге-Кутта 4-го порядкаtl012345678910У|h=1h=0.51.00000.64581.88804.97881.00000.63291.86594.951710.012917.025726.030537.032350.032965.033282.03339.983316.994925.999237.000750.001365.001582.0016h=0.1ИстинноерешениеОшибкаh=1h=0.5h=0.11.00000.63211.86474.95029.981716.993325.997536.999149.999764.999982.00001.00000.63211.86474.95029.981716.993325.997536.999149.999764.999982.00000.00000.01370.02330.02860.03120.03240.03300.03320.03320.03330.03330.00000.00080.00120.00150.00160.00160.00170.00160.0016L 0.00160.00160.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00004 6 8 • Глава 14.
Дифференциальные уравненияКак видно из таблицы, при длине шага, равной 0.1 (используемой системой по умолчанию), решения оказываются точными до четвертого знака. Уменьшение длины шага до0.01 увеличивает точность до седьмого знака. К сожалению, подобная закономерностьреализуется далеко не для всех функций. Выбор желаемого уровня точности остаетсяза пользователем, но зачастую оценить величину ошибки при той или иной длине шагапросто невозможно. В этом случае альтернативным оказывается метод Рунге-Куттас переменным шагом (адаптивный метод), лежащий в основе работы функций rkadaptи Rkadapt. В каждой узловой точке интервала этот алгоритм находит два варианта решения, используя схемы Рунге-Кутта 4-го и 5-го порядков.
Если разность между полученными значениями оказывается меньше допустимой величины ошибки (которуюможно заранее задать через системную переменную TOL), то длина шага увеличивается. Если же разность превысит ошибку, то следующую точку интервала алгоритмопределит с меньшей длиной шага. В случае равенства обоих параметров шаг остается постоянным. Адаптивный метод целесообразно использовать для нахождения неоднородных решений, меняющихся с различной скоростью на том или ином участке интервала. Это позволит сократить время расчета, поскольку для некоторых функцийхороший результат может быть получен и при использовании алгоритмом небольшого количества шагов.Метод Булирша-ШтераАлгоритм Булирша-Штера (Bulirsch-Stoer), реализованный во встроенных функциях bulstoer и Bulstoer, дает хорошие результаты для гладких функций решения. В егооснове лежит численное приближение решения на промежутке [х;х+Н] методом рациональной экстраполяции.
Серия значений, полученных на интервале интегрирования,экстраполируется в конечной точке интервала к истинному решению полиномом, приближающим функцию (рис. 14.20).хх+НРис. 14.20. Иллюстрация метода Булирша-ШтераПриближение оказывается тем точнее, чем больше шагов было выбрано внутри интервала. Понятно, что при стремлении их количества к бесконечности оно превращаетсяв точное аналитическое решение, но при этом неограниченно возрастает и время расчета. Сократить его можно, задав требуемую точность решения TOL через параметр асе.14.2. Численное решение ОДУ в форме задачи Коши* 469На промежутке последовательно выбираются п=2, 4, 6,8,10 и т.
д. шагов интегрирования с постоянно уменьшающейся длиной до тех пор, пока разность между полученными на конце интервала значениями не окажется равной TOL. Последний результат заносится в матрицу решений, выбирается следующий промежуток, и все операцииповторяются заново. Проанализируем конкретный пример, чтобы оценить, каким образом влияют на точность решения выбранное алгоритмом количество шагов и величина TOL.Пример 14.19. Влияние количества шагов и величины TOL на точностьалгоритма Булирша-ШтераГi ^Sys(t,y):=уО:=^-yo-15yj3x:=bulstoer(yO, 0,10,l(f ,Sys, 10,0.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
