Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 102
Текст из файла (страница 102)
Чтобы существовала возможность оценки точности этих функций в зависимости от различных параметров,расчет будем производить для ОДУ второго порядка, решение которого элементарнонаходится аналитически. Кроме того, в нем мы попытаемся продемонстрировать всеизложенные выше трудности при применении рассматриваемых функций.Пример 14.16. Функции rkadapt и bulstoerРешаемое уравнение второго порядка с начальными условиями:у"(х)+ 10у'(х)-х+ 10у(х) = 0У'(0) = 0у(0) = 1Его точное аналитическое решение:.
,-5ху(х) := еПреобразование ОДУ второго порядка в систему ОДУ в векторной форме (переносим в правуючасть уравнения все слагаемые, кроме содержащего вторую производную, и вводим заменыУо=У(х)> У'о"У,):УD(x,y) :=У0:=Численное решение полученной системы:rl :=rkadapt^O,O,l,Ю"4,D, 100,0.000l)bl :=bulstoer(yO,0,l,10~4,D, 100,0.000l)Оценка количества шагов, потребовавшихся для поиска решения с заданной точностью:(l)last( r l ) = 8(l)last(bl ) = .4 5 8 • Глава 14. Дифференциальные уравненияМатрицы решений:*О1ОО.О11-0.10.060.982-0.5890.310.618-1.917rl = 0.5030.7030.282-1.4180.084-0.5930.8290.032-0.2670.9290.013Ы =О1ОО.О11-0.10.110.941-1.0350.6150.151-0.9271-36.696x10-0.067-0.124316.7ЦхЮ"-0.067,чОценка погрешности исходя из максимального отклонения от аналитического решения:max(|rl ( 1 ) - у ( г 1 < 0 > ) | ) = L855 х 10" 4m a x ( | b l ^ _ у ( Ь 1 ® ) | ) = 4.365 х Ю"5Оценка влияния на точность решения уменьшения величины асе на три порядка:r2 :=rkadapt(yO,0,1,10~ 7 , D , 100,0.000l)last^i-2Ь2 :=bulstoer(yO,0,1,10" 7 , D , 100,0.000l)) = 26Iast(b2 ) = 57/JOhтах^^-у^^З^ЗЗхЮ--y(b2-9) | ) = 8.803x10Оценка значимости правильного согласования параметров:r3:=rkadapt(yO,O,l,lO- 15 ,D,lOO,o)Iast(r3^)=99^г398,0=0 Л 7 8- y(r3 < 0 > )|) = 2.109 х 10- ИЬЗ :=bulstoei/yO,0,l, \6~ 1 5 , D , 100,о)1 8 8 ^ 3 ^ ) = 11mЬ 3 1 0 Q = 0.957a x ( | b 3 0 > -y ( b 3 < 0 > ) | ) = 3.015 х 10.-14Решения системы ОДУ, полученные с различной точностью с помощью функцийrkadapt и bulstoer, представлены на рис.
14.12.Рис. 14.12. Решения системы ОДУ, полученные с различной точностью с помощью функцийrkadapt и bulstoer14.2. Численное решение ОДУ в форме задачи Коши * 4 5 9Изучая приведенный пример, обратите внимание на следующее.• Точность метода Булирша-Штера в случае нашей системы оказалась выше, чему адаптивного алгоритма Рунге-Кутта при том, что и шагов он использовал меньше. Этот факт можно объяснить тем, что искомые решения являются гладкими и неосциллирующими функциями.
При решении же более сложных систем эффективность адаптивного алгоритма может оказаться выше, чем метода Булирша-Штера.•Чтобы найти решение при всех рассмотренных уровнях точности, понадобилосьнеожиданно мало шагов (особенно это касается функции bulstoer). Учитывая это,не спешите делать количество разбиений при решении ОДУ с помощью стандартных функций большим — для получения результата обычной точности (0.001) дляэтого бывает вполне достаточно всего несколько десятков шагов.•Если вы хотите построить график решений системы ОДУ, сделайте величину асе1015поменьше (10~ -10~ ).
В этом случае система сгенерирует достаточно точек длязадания гладкой кривой. При обычной же точности по полученным данным вы сможете построить лишь весьма грубую ломаную.•Очень важно правильно согласовать параметры точности, максимального количествашагов и минимальной их длины. Так, в приведенном примере при асс=10~15 и отсутствии ограничений на размер шага адаптивный алгоритм смог просчитать функциирешения только на 1/10 длины интервала (см.
правый график рисунка из примера 14.16). Для того же, чтобы получить с помощью него корректное решение, лимитколичества шагов должен быть увеличен до 2000. Зато при том же наборе параметров функция bulstoer нашла результат всего при 11 разбиениях — более чем в 100 разменьшем их количестве, чем понадобилось rkadapt. Этот факт как нельзя лучше демонстрирует, как сильно определяет успех расчетов правильный выбор алгоритма.Чтобы проиллюстрировать, как различается шаг на промежутке при адаптивном егоопределении, построим, исходя из полученных в примере 14.16 значений, диаграмму(рис.
14.13).i := 0.. Iast(r2 <i;> )j := 0.. last(j>2^)110.8онИ---" •••••l0.5Рис. 14.13. Переменный шаг, использованный функциями rkadapt и bulstoerИз приведенной диаграммы следует, что шаг, используемый численными методами,зависит как от скорости изменения функций системы, так и от их величин.Чтобы продемонстрировать все возможности Mathcad в решении систем ОДУ, рассмотрим на примере нахождения траектории ракеты, насколько эффективно справляется с поставленной задачей вычислительный блок Given-Odesolve.4 6 0 * Глава 14.
Дифференциальные уравненияПример 14.17. Движение ракеты в поле тяготения небесных телСогласно открытому Ньютоном закону всемирного тяготения, все тела притягиваются другк другу с силой, прямо пропорциональной произведению их масс и обратно пропорциональнойквадрату расстояния между ними (в данном случае объекты рассматриваются как материальныеточки, то есть размерами тел по сравнению с расстоянием между ними можно пренебречь):В нашем случае m — масса ракеты, М — масса планеты, R — расстояние между ними, G — гравитационная постоянная (чтобы не оперировать большими величинами масс и расстояний, не будем учитывать ее в примере, поскольку на характер зависимости это никак не повлияет). С другой стороны, второй закон Ньютона гласит, что сила есть не что иное, как произведение массытела на его ускорение:F = maПриравняв выражения и отбросив G, получим:МИспользуя теорему Пифагора, несложно показать, что расстояние R между двумя точками с координатами (х, у, z) и (xl, yl, zl) определяется следующим выражением:l-x)2 + (yl-y) 2 +(zl-z) 2Для упрощения расчетов найдем проекции вектора ускорения а на координатные оси х, у, г.
Чтобы наглядно представить все рассуждения, начертим вспомогательную схему (рис. 14.14).Рис. 14.14. Схема нахождения проекции вектора ускорения а на ось хПроекция расстояния R между ракетой и планетой М является катетом прямоугольного треугольника Мх1х (см. рис. 14.14), а само расстояние — гипотенузой. Поскольку косинус угла естьотношение прилежащего катета к гипотенузе, то проекцию вектора а мы можем найти какax=a-cos(a)Подставим в полученную формулу выражение для а и воспользуемся определением косинуса,в результате получимМ"3"R(xl-x)M-(xl-x)2:1-х)+(yl-y)+(2l-z) 2 ]Абсолютно аналогичные выражения можно записать и для проекций вектора а на оси у и z.Полученная нами зависимость описывает траекторию полета ракеты в поле тяготения одной:планеты.
В случае же нескольких планет необходимо учитывать принцип суперпозиции грави-14.2. Численное решение ОДУ в форме задачи Коши • 4 6 1тационных полей. Согласно этому принципу гравитационное поле, создаваемое массой, действуетвне зависимости от наличия других масс. Другими словами, гравитационные поля, возбуждаемые несколькими небесными телами, накладываются друг на друга, оставаясь при этом неизменными, и управляют движением как искусственных, так и естественных тел в космическом пространстве. В нашем примере поле тяготения каждой планеты будет вносить вклад в ускорениеракеты. Сумма же этих вкладов и будет результирующим ускорением.
В дифференциальнойформе для составляющей вектора ускорения по оси х сказанное выше запишется как:M r (xl-x)ipdtM2-(x2-x)+IT5 +Mn-(xn-x)+'"R"1По такому же принципу составляются уравнения и для проекций вектора а по осям у и z.Теперь, когда даны необходимые теоретические пояснения, приступим к решению полученнойсистемы дифференциальных уравнений.Траектория полета определяется рядом факторов, важнейшие из которых — начальная скоростьракеты, масса, координаты планет, а также их количество. Варьируя эти параметры, рассмотримособенности траектории ракеты в наиболее интересных случаях.Начнем рассмотрение задачи с простейшего случая, когда ракета, имеющая нулевую начальнуюскорость, движется в направлении одной планеты.Если начальная скорость ракеты равна нулю, то она будет двигаться в направлении планетыпо прямой линии.
Аналогичная траектория будет наблюдаться и в том случае, когда начальнаяскорость ракеты направлена точно к центру планеты.Изучая траекторию движения ракеты, мы будем увеличивать количество планет до пяти, поэтому для задания масс и координат планет по мере их увеличения воспользуемся таблицами вводаданных (в общем случае количество планет может быть произвольным).В данном примере приводятся дифференциальные уравнения для универсальной модели, описывающей движение ракеты в поле тяготения п планет. Если же задать в таблице параметры одной планеты, то решение найдено не будет, поскольку для корректной работы модели массыи координаты должны быть представлены в форме векторов, что возможно, когда минимальноеколичество планет равно двум.
Поэтому добавим в систему еще одну планету малой массы, влиянием которой на движение ракеты можно пренебречь.М:=х:=z:=_у:=00301-Ю11-10-87020010020500101Определяем координаты исходной точки и составляющие скорости ракеты по осям х, у, z.Vx:=0Vy:=0x0:=0y0:=0Vz:=0z0:=0Вводим ключевое слово Given, начальные условия и общие уравнения движения.Givenх(0) = х0x'(0) = VxУ(О) = уОy'(0) = Vyz(0) = z0z'(0) = Vz1030462 *Глава 14. Дифференциальные уравненияlast(x)2— 2 x(t) =dtlast(y)M,(y.- y(t))i=0,2last(z)—2i=0M.YZ.r.пЗНаходим решение системы,t, 150,4000R:=OdesolveЧтобы воспользоваться функцией CreateSpace для визуализации решения, зададим входящийв нее параметр — вектор-функцию F.: RХ:== oY:=RF(t) :=Z:=: = R2Y(t)В том, что система была решена нами верно, легко убедиться, взглянув на рис. 14.15.(x,y,z),CreateSpace(F,O,:5O,200Cf)Рис.
14.15. Траектория движения ракеты в гравитационном поле одной планеты при нулевойначальной скорости14.2. Численное решение ОДУ в форме задачи Коши * 4 6 3Рассмотрим теперь движение ракеты, когда ее начальная скорость не направлена на планету.Когда начальная скорость ракеты не направлена на планету, траектория полета искривляетсяпритяжением. Если скорость ракеты не превышает некоторой величины, то по закону Кеплератраектория будет представлять собой эллипс, который лежит в плоскости, проведенной черезначальное направление скорости и центр планеты. Центр же притяжения будет находиться в одном из фокусов эллипса. Разумеется, чем больше начальная скорость, тем больше будет ось орбиты.В случае слишком большой скорости ракета сможет уйти от центра притяжения.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
