Численный пример. В качестве иллюстрации рассмотрим примене-
ние некоторых из описанных выше методов для решения простейшей фи-
зической задачи о движении частицы массой mпод действие силы упруго-
сти F=-kx. Определим зависимость положения частицы от времени, если в
начальный момент времениt
0
=0 ее положение и скорость имеют заданные
значенияx(0)=x
0
,
(0)=
0
. Движение частицы описывается уравнением
Ньютона:
(18)
Это обыкновенное дифференциальное уравнение второго порядка,
которое можно заменить эквивалентной системой двух уравнений первого
порядка:
,
(19)
где v – новая функция, имеющая в данном случае смысл скорости.
14
Уравнения (19) необходимо дополнить начальными условиями (постанов-
ка задачи Коши):x(0)=x
0
, υ(0)= υ
0
.
Решение этой задачи хорошо известно: частица будет совершать
гармонические колебания относительно положения равновесия (x=0) по
законуx(t)=Asin(ωt+φ
0
), где A – амплитуда колебаний,
– частота колеба-
ний, равная
, φ
0
– начальная фаза колебаний. Наличие точного анали-
тического решения этой задачи дает нам возможность оценить точность
приближенных численных результатов.
Перейдем к обсуждению численного решения системы уравнений
(19). Предположим, что решение уравнений (19) ищется на интервале
времениtϵ[0,T]. Введем на этом интервале дискретную сетку узлов с
шагомΔt: {t
i
=i·Δt, Δt=T/N, i=0, 1, …, N}
Применение метода Эйлера (8) для решения системы уравнений (19)
приводит к простым расчетным формулам:
x
i+1
=x
i
+ υ
i
Δt, υ
i +1
= υ
i
-k/m · x
i
·Δt, i=0, 1, …, N(20)
Применение формулы модифицированного метода Эйлера (10) сво-
дится к последовательному вычислению значений функций x и v сначала с
шагом 1/2, с целью определения наклона интегральных кривых в середине
отрезка интегрирования, а затем с шагомΔt:
,
,
Исключая из последних формул промежуточные значенияx
i+1/2,
, получим:
,
, (21)
i=0,1,…,N
Можно показать, что применение метода Эйлера с пересчетом (12) к
данной системе уравнений дает результат, в точности совпадающий с (21).
И, наконец, применим для решения системы уравнений (19) формулы ме-
тода Рунге-Кутта четвертого порядка (14):
ℎ
ℎ
ℎ
ℎ
ℎ
ℎ
ℎ
Ниже в таблице 1.1 приведены результаты численных расчетов при
начальных условияхx(0)=1,v(0)=0 для различных моментов времени.Масса
частицы m и коэффициент упругости k приняты равными единице. Шаг
интегрированияΔt во всех методах равен 0,1.
15
Таблица 1.1
Результаты численных расчетов
Момент
времени
t
Точное
значение
x(t)
Метод Эй-
лера
Модифицированный
метод Эйлера
Метод
Рунге-Кутта
IV
0
1
1
1
1
2
-0,41615
-0,45302
-0,41927
-0,41615
4
-0,65364
-0,80974
-0,64892
-0,65365
6
0,96017
1,28642
0,96363
0,96017
8
-0,14550
-0,17751
-0,15880
-0,14549
10
-0,83907
-1,40885
-0,83095
-0,83908
Из приведенных результатов видно, что метод Эйлера обладает наи-
большей погрешностью, которая, к тому же, со временем сильно растет,
т.е. характеризуется склонностью к накапливанию. Модифицированный
метод Эйлера дает существенно лучший результат. На всем
интервалеtϵ[0,10].Относительная погрешность его результатов не превос-
ходит 10%. Метод Рунге-Кутта показывает очень хорошую точность своих
результатов с максимальной относительной погрешность не более 0,01%.
Выбор шага интегрирования и контроль за точностью вычислений.
До сих пор мы не обсуждали вопрос выбора шага интегрирования, а ведь
именно от величины шага зависит точность получаемого решения и время,
затрачиваемое на его получение. Из приведенных ранее рассуждений о по-
рядке точности методов решения ОДУ можно сделать общий вывод о том,
что для повышения точности следует брать меньший шаг. Однако на прак-
тике ситуация не столь однозначна. Уменьшение шага интегрирования
приводит к увеличению времени вычислений. И что более важно, слишком
малые значения шага могут привести не к повышению точности, а, наобо-
рот, к увеличению погрешности в силу накапливания вычислительной
ошибки. В тоже время, выбор слишком большого шага интегрирования
может привести не только к большой погрешности, но и к получению аб-
солютно неверного результата. Поэтому выбор шага это всегда определен-
ный компромисс между точностью и временем.
Таким образом, мы обозначили проблему выбора такого значения
шага интегрирования, при котором бы обеспечивалась требуемая точность
вычисления и умеренные затраты машинного времени. И здесь, прежде
всего, нужно определить критерий, по которому можно судить о точности
получаемых результатов. В рассмотренном выше численном примере та-
ким критерием являлось сравнение приближенного результата с точным.
Однако, если точный результат неизвестен, что и бывает в подавляющем
большинстве случаев, то применение такого критерия оказывается невоз-
16
можным.
Одним из возможных критериев точности может служить сравнение
приближенных результатов в каждом узле, полученных при разных шагах
интегрирования, например h и h/2. Если величина
сравнима с заданной погрешностью вычислений, то шаг можно увеличить;
в противном случае, когда указанная величина слишком велика, значение
шага следует уменьшить. Используя эту оценку, можно построить методы
с автоматическим выбором шага и контролем точности на протяжении все-
го времени вычислений. Такие алгоритмы называют адаптивными, т.е.
подстраивающимися под условия конкретной задачи.
Другой подход опирается на использование особенностей решаемой
задачи. В частности, при решении многих физических задач может суще-
ствовать такой параметр, например, полная энергия, который должен со-
хранять свое значение. В ходе вычислений, постоянно отслеживают значе-
ние этого параметра и, в случае его изменения, выходящего за пределы до-
пустимой погрешности, корректируют значение шага интегрирования.
Что касается выбора начального (пробного) значения шага, то здесь,
к сожалению, не существует универсального рецепта. И в каждом кон-
кретном случае шаг выбирается исходя из характерных параметров решае-
мой задачи.
Do'stlaringiz bilan baham: |