Экспериментальное исследование свойств методов Рунге-Кутты
МИНИСТЕРСТВО НАУКИ И ОБРАЗОВАНИЯ УКРАИНЫ
Севастопольский национальный технический университет
Кафедра технической кибернетики
КУРСОВАЯ РАБОТА ПО ДИСЦИПЛИНЕ
«Вычислительные методы» на тему:
«Экспериментальное исследование свойств методов Рунге-Кутты»
Выполнила: студентка гр. А-31д
Воротилова Я.М.
Проверил: Мирянов В.И.
Севастополь
2004
СОДЕРЖАНИЕ
ВВЕДЕНИЕ
1 ПОСТАНОВКА ЗАДАЧИ
1.1 Приведение к нормальной форме Коши
1.2 Метод Рунге-Кутты второго порядка
2 ОПИСАНИЕ ПРОГРАММНЫХ МОДУЛЕЙ
2.1 Основная программа
2.2 Функция вычисления точного решения
2.3 Процедура вычисления правых частей системы уравнений в нормальной форме Коши
2.4 Процедура RK2
2.5 Процедура RK4
3 ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ МЕТОДОВ РУНГЕ-КУТТЫ
3.1 Анализ влияния величины шага на точность интегрирования методами Рунге-Кутты второго и четвертого порядка
3.2 Проверка гипотезы Рунге
3.3 Исследование поведение ошибки интегрирования как функции независимой переменной для обоих методов Рунге-Кутты при различных значениях шага
3.4 Сравнительный анализ эффективности методов Рунге-Кутты при различных требованиях к точности вычисления
ЗАКЛЮЧЕНИЕ
БИБЛИОГРАФИЯ
ПРИЛОЖЕНИЕ А
ПРИЛОЖЕНИЕ Б
ПРИЛОЖЕНИЕ В
ВВЕДЕНИЕ
Настоящая курсовая работа посвящена опытному исследованию свойств методов Рунге-Кутты и реализации на персональных компьютерах численных методов приближенного интегрирования ОДУ, наиболее часто применяющихся в практике моделирования и проектирования СА и У. Экспериментальные исследования проводятся с помощью составленных и отлаженных программ интегрирования обыкновенных дифференциальных уравнений на ЭВМ.
Задание предполагает:
закрепление теоретических навыков и знаний в вопросе о проблематике интегрирования ОДУ и численного решения задачи Коши методом Рунге-Кутты, изучение их основных свойств (точность, эффективность, устойчивость) и основных характеристик данных свойств (локальная и глобальная алгоритмические ошибки, порядок метода, ошибка вычисления и т.п.) ;
приобретение основных навыков составления и отладки процедур и функций интегрирования на основе методов Рунге-Кутты и программ интегрирования систем дифференциальных уравнений с использованием все тех же процедур и функций;
проведение опытных исследований зависимости точности, эффективности и устойчивости алгоритмов интегрирования от величины шага интегрирования и порядка метода Рунге-Кутты на ЭВМ.
В различных сферах технических и даже экономических отраслей приходится достаточно часто сталкиваться с математическими задачами, для которых не представляется возможным описать точное решение классическими методами или сие решение выражено крайне неудобочитаемыми соотношениями, которые представляют из себя неприемлемую для мозга пищу, не говоря уже об использовании или реализации на практике.
Разрабатываемые вычислительной математикой численные методы носят в основном ориентировочный характер, однако они позволяют получить итоговый числовой результат со сносной для практических нужд точностью. Численные методы представляют собой алгоритмы вычисления приблизительных значений искомого решения на определенной сетке значений аргумента. При определенных условиях значения аргумента могут являться точными.
Численные методы не позволяют найти общее решение: полученное решение является частным. Но одним из многочисленных плюсов данных методов можно назвать высокую степень применимости к обширным классам уравнений и всем типам вопросов и заданий к ним. Посему с появлением электронных вычислительных машин численные методы стали одними из основных технологий решения определенных практических задач решения ОДУ.
Большую значимость имеет вопрос о верности вычислений на ЭВМ, поскольку при практической реализации имеет место обширный объем обрабатываемой подсчитываемой информации и погрешности могут достаточно сильно исковеркать конечный результат, принимаемый нами за действительный с «поправками на ветер». Кроме сказанного оценка точности численного метода немаловажна и потому, что увеличить точность в некоторых пределах можно за счет увеличения объемов вычислений, а уменьшить временные затраты при решении задачи - за счет снижения точности получаемого результата.
Для
понижения погрешности методов
интегрирования ОДУ, использующего
разложения искомого решения в ряд
Тейлора, необходимо принимать во внимание
большее количество членов ряда. При
всем при этом появляется потребность
аппроксимации производных правых частей
ОДУ. Ключевая идея методов Рунге-Кутты
заключается в том, что производные
аппроксимируются через значения функции
в точках на интервале
,
которые выбираются из условия наибольшей
близости алгоритма к ряду Тейлора. В
зависимости от старшей степени
,
с коей учитываются члены ряда, построены
всевозможные вычислительные схемы
Рунге-Кутты разных порядков точности.
Среди достоинств схем Рунге-Кутты не следует обходить во внимании:
удобоваримую точность;
одноступенчатость,
то есть дабы найти
,
необходима информация лишь о предыдущей
точке
;
координирование
с рядом Тейлора вплоть до членов порядка
,
где степень
неодинакова для различных методов и
именуется порядком метода;
отсутствие
необходимости вычисления производных
от
,
причем накладывается требование
вычисления всего-навсего самой функции.
Собственно
благодаря вышеуказанному свойству c)
методы Рунге-Кутты предпочтительней
рядов Тейлора для реализации на практике.
Тем не менее поводов для веселья мало,
ибо перед нами стоит нелегкая задача
неоднократного вычисления функции
при неодинаковых значениях
и
для вычисления последующей точки
решения. Это Богом дарованное наказание
за преподнесенную нам численным методом
поблажку, заключающуюся в отсутствии
какой бы то ни было надобности вычисления
иной раз весьма громоздких производных,
но трудностей боятся кто угодно, только
не мы.
1 ПОСТАНОВКА ЗАДАЧИ
1.1 Приведение к нормальной форме Коши
Нормальной формой Коши принято
называть общую форму записи ОДУ, то есть
представление в виде системы
дифференциальных уравнений первого
порядка:
(1)
ДУ второго порядка, заданное согласно варианту №3 имеет вид:
(2)
Задание
предполагает нахождение решения
на интервале
при следующих начальных условиях:
(3)
Для решения ДУ его просто необходимо представить согласно нормальной формы Коши. Для этого руководствуемся следующими обозначениями:
(4)
В итоге имеется система ДУ первого порядка вида:
(5)
Произведя все вышеописанные манипуляции над заданным в варианте уравнением, получим следующую систему:
(6)
Система (6) есть решение уравнения (2).
1.2 Метод Рунге-Кутты второго порядка
В методах Рунге-Кутты интеграл заменяется линейной комбинацией значений подынтегральной функции, вычисленных при разных значениях аргумента:
(7)
Метод Рунге-Кутты представим в виде:
(8)
Из вышеуказанных общих формул (8) получают формулы метода Рунге-Кутты 2-ого порядка m=2;
(9)
Для определения метода необходимо
найти значения вещественных коэффициентов:
.
Для этого интеграл, заменяемый
линейной комбинацией значений
подынтегральной функции, вычисленных
при разных значениях аргумента, можно
представить как:
(10)
А его, в свою очередь, можно представить рядом Тейлора:
(11)
где
-
сумма элементов ряда Тейлора, степень
которых не ниже 3.
Осталось найти неизвестные
значения
(12)
В результате таких бесхитростных манипуляций получаем искомый ряд Тейлора:
(13)
Приравняем коэффициенты при
одинаковых степенях
в выражениях
(11) и (13). В итоге получим систему уравнений вида:
(14)
Из свойств системы (14) следует
отметить, что она не обладает единственным
решением. При
значение
,
значение
,
а
(15)
Подставив полученные коэффициенты в соотношение (8), получаем следующие формулы метода Рунге-Кутты 2-ого порядка:
(16)
2 ОПИСАНИЕ ПРОГРАММНЫХ МОДУЛЕЙ
Составленная в ходе курсовой работы программа вычисляет решения дифференциального уравнения, с предварительно заданными начальными условиями. Интегрирование происходит согласно двум методам: Рунге-Кутты второго и четвертого порядков.
Программа состоит из следующих модулей:
Основная программа;
Процедура вычисления точного решения ДУ;
Процедура вычисления правых частей;
Процедура выполняющая шаг интегрирования методом Рунге-Кутты 2-ого порядка;
Процедура выполняющая шаг интегрирования методом Рунге-Кутты 4-ого порядка.
2.1 Основная программа
Блок программы осуществляет следующие операции:
запрашивает у нерадивого пользователя величину шага интегрирования и шаг вывода на экран;
вычисляет количество шагов;
с заданным шагом вызывает процедуры интегрирования методом Рунге-Кутты 2-ого и 4-ого порядков на отрезке интегрирования;
вычисляет погрешность и оценку погрешности интегрирования;
выводит замечательные результаты работы программы с заданным шагом вывода на экран.
Для простоты понимания укажем следующие переменные, содержащиеся в программе:
h - шаг интегрирования. Вводится нерадивым пользователем с клавиатуры;
n – число шагов интегрирования;
h_screen - шаг вывода результатов на экран. Вводится нерадивым пользователем с клавиатуры;
i_screen – счётчик вывода результатов на экран. Когда i_screen > h_screen, то происходит вывод результатов и обнуление i_screen;
i, j – переменные, используемые циклом;
e2, e4– ошибки интегрирования для методов Рунге-Кутты 2-ого и 4-ого порядков соответственно. Подсчитываются из соотношения(1):
(1)
e2max, e4max – оценки погрешностей интегрирования для методов Рунге-Кутты 2-ого и 4-ого порядков соответственно. Подсчитываются из соотношения(2):
(2)
t – значения независимой переменной;
t0, tf – пределы интегрирования
y2, y4 – вектора решения для методов Рунге-Кутты 2-ого и 4-ого порядка соответственно в узле t>k>> >;
outfile– переменная файлового типа. Определена для вывода результатов в текстовой файл;
name – переменная строкового типа. Используется для передачи имени файла.
Текст основной программы приведён в приложении А, схема в приложении Б.
2.2 Функция вычисления точного решения
function clearsolve (t: real): real
Функция предназначена для вычисления точного решения для дифференциального уравнения по формуле (3):
(3)
Текст функции приведен в приложении 2, схема в приложение 7.
2.3 Процедура вычисления правых частей системы уравнений в нормальной форме Коши
procedure right(t: real; var x,f: vector_n);
Процедура вычисляет правые части системы однородных дифференциальных уравнений в нормальной форме Коши по формуле (4):
(4)
Текст процедуры приведен в приложении А, а схема в приложение Б.
2.4 Процедура RK2
procedure RK2(t: real; h: real; var x: vector_n);
Укажем формальные параметры:
t – независимая переменная ;
h – шаг интегрирования;
x – массив решений. При входе в процедуру решение в текущем узле интегрирования, при выходе в следующем.
Процедура RK2 выполняет шаг интегрирования системы ОДУ методом Рунге-Кутты 2-ого порядка из соотношения (5):
(5)
где
(6)
Процедура обращается к процедуре
вычисления правых частей right
с различными параметрами для вычисления
и
(6).
Затем с Божьей помощью (5) считает значение
.
Текст процедуры приведен в приложении А, схема в приложение Б.
2.5 Процедура RK4
procedure RK4(t: real; h: real; var x_4: vector_n);
Формальные параметры:
t – независимая переменная ;
h – шаг интегрирования;
x – массив решений. При входе в процедуру решение в текущем узле интегрирования, при выходе в следующем.
Процедура RK4 выполняет шаг интегрирования системы обыкновенных дифференциальных уравнений (1.1.2) методом Рунге-Кутты 4-ого порядка (7).
(7)
где
(8)
Процедура четыре раза обращается
к процедуре вычисления правых частей
right с разными параметрами
для вычисления
,,
,
(8).
Затем с Божьей помощью (7) считает значение
.
Текст процедуры приведен в приложении А, схема в приложение Б.
3 ЭКСПЕРИМЕНТАЛЬНОЕ ИССЛЕДОВАНИЕ МЕТОДОВ РУНГЕ-КУТТЫ
3.1 Анализ влияния величины шага на точность интегрирования методами Рунге-Кутты второго и четвертого порядков
Плоды деятельности на ПЭВМ приводятся в приложении В. Результат упорного труда программы представлены в графическом виде (рисунок 1).
Рисунок 1 – Зависимость оценки e2 от шага интегрирования
Рисунок 2 – Зависимость оценки e4 от шага интегрирования
На рисунках изображены зависимости оценки погрешности интегрирования от величины шага интегрирования для обоих методов. Из него видно, что практические результаты соответствуют теоретическим положениям, но не совсем.
При шаге от 0.1 до 1 бесспорно влияние погрешности интегрирования как в методе второго, так и в методе четвертого порядков. Ошибка на данном интервале начинает лавинообразно возрастать, что связано с нарушением устойчивости алгоритма.
С дальнейшим уменьшением шага до 0.1 – 0.001 величина погрешности уменьшается за компанию, и наблюдается достаточно большая точность вычислений.
Дальнейшее уменьшение шага (менее 0.001) вызывает увеличение полной ошибки, а также скверное её поведение. Это связано с возрастанием влияния ошибки вычислений из-за увеличения количества вычислений, необходимых для получения решения. На фоне уменьшения алгоритмической погрешности решающую роль играет погрешность вычислений, которая представляет собой сумму всех ошибок округления при реализации данного метода на конкретной ПЭВМ.
Из двух методов 2-го и 4-го порядков при одинаковых значениях шага точнее метод четвертого порядка, но при уменьшении шага, точность методов постепенно выравнивается.
3.2 Проверка гипотезы Рунге
Согласно гипотезе, суммарная погрешность алгоритма при интегрировании ДУ с постоянным шагом пропорциональна величине шага в степени, равной порядку метода.
(1)
Результаты вычислений на ПЭВМ приводятся в приложении В. Результат работы программы наглядно представлены в графическом виде (рисунок 2).
Рисунок 3 – Зависимость отношения оценки погрешности к величине шага интегрирования в степени, равной порядку метода от шага интегрирования
Гипотеза Рунге экспериментально подтверждается только при таких значениях шага (1-0.01), где погрешность вычисления ПЭВМ влияет на результат в благоразумных пределах. При дальнейшем уменьшении шага пропорциональность скрывается из виду. Это можно истолковать тем, что предположение Рунге не учитывает влияния ошибки вычисления ПЭВМ на полученный результат.
3.3 Исследование поведение ошибки интегрирования как функции независимой переменной для обоих методов Рунге-Кутты при различных значениях шага
Для проведения своих изысканий нам необходимо выбрать 3 значения шага в начале, середине и конце диапазона изменения шага и построить графические зависимости решений от независимой переменной и погрешностей интегрирования от независимой переменной.
Результаты вычислений на ЭВМ приводятся в приложении В. Результат работы программы наглядно представлены в графическом виде (рисунок 3 и 4).
Рисунок 4 – График точного и приближенных решений (h=0.001, 0.1, 0.5) дифференциального уравнения
На рисунке 3 слегка, но все-таки заметны отклонения приближённых решений при разной величине шага от точного, что подтверждает высокую точность получаемого приближенного решения относительно точного.
Рисунок 5 – График зависимости величины ошибки интегрирования от независимой переменной при h=0.001, 0.1, 0.5
При шаге, численно равном 0.5 из рисунка 5 очевидно, что различия между погрешностями методов 2-ого (E5) и 4-ого (E6) порядков довольно значительны, что объясняется более высокой точностью метода 4-ого порядка.
Разница ошибок интегрирования для обоих методов при шаге интегрирования 0.1 остается достаточно высокой, но следует подметить, что точность метода Рунге-Кутты 4-ого порядка при данном значении шага выше, чем при меньшем h=0.001 (это обусловлено тем, что алгоритмическая ошибка стремится к нулю, а вычислительная ещё не показывает свой отвратительный характер).
При величине шага интегрирования h=0.001 ошибки интегрирования (E1 и E2), равно как и при больших h возрастают с ростом независимой переменной, что слава Богу не противоречит теории, однако прослеживается скачкообразный характер графиков, что связано с внушительной лептой вычислительной ошибки в общую погрешность.
3.4 Сравнительный анализ эффективности методов Рунге-Кутты при различных требованиях к точности вычисления
Для анализа надлежит оценить затраты машинного времени на интегрирование ДУ в зависимости от величины погрешности интегрирования.
Так как для довольно непростых уравнений львиная доля времени уходит на вычисление правых частей уравнений, то в качестве оценки затрат машинного времени можно принять количество вычислений правых частей при нахождении решения на всём отрезке интегрирования, что есть произведение количества шагов интегрирования на порядок метода. Зависимость оценки ошибки интегрирования от количества шагов интегрирования определена.
Замечательные результаты вычислений на ПЭВМ приведены все в том же приложении В. Результат работы программы представлены графически (рисунок 5).
Рисунок 6 – График зависимости величины ошибки интегрирования от количества вычислений правой части
Из близлежащего графика практически определяется метод для различных требованиях к точности и времени работы нашей замечательной программы.
График наглядно демонстрирует, что наличествует такой предел шага, ниже которого программа гоняет свои байты зря, так как оценка погрешности растёт из-за ошибки вычисления ПЭВМ. Стало быть, на практике следует выбирать определенный промежуток шага, желательно в котором алгоритм устойчив.
ЗАКЛЮЧЕНИЕ
При реализации на практике задания для курсовой работы, заключающегося в интегрирования ОДУ, была составлена и отлажена программа, приведенная в приложении А. С помощью данной программы проведена серия опытных исследований свойств методов Рунге-Кутты второго и четвёртого порядков.
При задании определенного
интервала значений шага интегрирования
ошибка интегрирования уменьшается с
уменьшением шага. Подтверждение сего
факта можно без труда найти в теории.
Однако нельзя обойти вниманием тот
труднооспоримый аргумент, что для
довольно незначительных значений
с уменьшением шага интегрирования
ошибке свойственно увеличивается. Это
связано с лавинообразным ростом числа
требуемых для получения решения
вычислений и c
увеличением ошибки вычислений.
На заданном интервале значений
(при условии, что на интервале ошибка
вычислений велика не до безобразия),
при одинаковых значениях шага
интегрирования метод Рунге-Кутты
четвёртого порядка имеет достаточно
малую ошибку вычисления относительно
ошибки метода Рунге-Кутты второго
порядка.
Следует отметить, что коэффициент
(формула (1) раздела 3) для таких значений
шага интегрирования, при которых можно
скрипя сердцем считать его неизменным,
имеет меньшее значение для метода
Рунге-Кутты четвёртого, нежели второго
порядка. Отсюда мораль: ошибка
интегрирования с помощью метода
Рунге-Кутты четвёртого порядка меньше
ошибки интегрирования при использовании
метода Рунге-Кутты второго порядка и
обусловлено это не только тем, что она
обратно пропорциональна величине шага
в четвёртой, а не во второй степени, но
и оттого, что коэффициент пропорциональности
при этом существенно меньше.
Метод Рунге-Кутты четвёртого порядка имеет значительно меньшую ошибку вычисления, нежели метод Рунге-Кутты второго порядка (при условии идентичных машинных затрат). При неизменной требуемой точности для метода Рунге-Кутты четвёртого порядка необходимы меньшие затраты машинного времени относительно метода Рунге-Кутты второго порядка, отчего при корректном выборе шага интегрирования метод Рунге-Кутты четвёртого порядка значительно результативнее, чем метод Рунге-Кутты второго порядка.
БИБЛИОГРАФИЯ
Самарский А.А., Гулин А.В. Численные методы .М.: Наука, ГРФМЛ, 1989.- 432с.
Бахвалов Н.С. Численные методы. М.: Наука, ГРФМЛ,1987.-600 с.
Ляшко И.И., Макаров В.Л. Интегральные методы вычислений. Киев, 1977, 408с.
Мантуров О.В. Курс высшей математики М.:В.Ш.-1991.-448с.
Маликов В.Т. Вычислительные методы и применение Киев: В.Ш.-1989.-213 с.
ПРИЛОЖЕНИЕ А
(Справочное)
Основная программа
program DU(input,output);
uses crt;
type vector_n = array [1..2] of real;
var t0, tf, k, j, n: integer;
t, yt, h, h_screen, e2, e4, e2max, e4max, i_screen: real;
y2, y4: vector_n;
name: string;
outfile: text;
begin
clrscr;
writeln('Please enter file name');
readln(name);
writeln('Please enter h,h_screen');
readln(h,h_screen);
clrscr;
writeln;
assign(outfile,name);
rewrite(outfile);
t0:=0;
tf:=10;
n:=round((tf-t0)/h);
y2[1]:=1;
y2[2]:=0;
y4:=y2;
e2:=0;
e4:=0;
e2max:=0;
e4max:=0;
t:=t0;
i_screen:=h_screen;
for k:=0 to n do
begin
yt:=clearsolve(t);
e2:=abs(yt-y2[1]);
e4:=abs(yt-y4[1]);
if e2>e2max then e2max:=e2;
if e4>e4max then e4max:=e4;
if i_screen>h_screen-0.00001 then
begin
yt:=clearsolve (t);
e2:=abs(yt-y2[1]);
e4:=abs(yt-y4[1]);
if e2>e2max then e2max:=e2;
if e4>e4max then e4max:=e4;
if i_screen>h_screen-0.00001 then
begin
writeln(' t=',t:6:3,'; yt=',yt:9:3,'; y2=',y2[1]:9:3,'; y4=',y4[1]:9:3,'; e2=',e2:6:3,'; e4=',e4:8:6);
writeln(outfile,' t=',t:10:6,'; yt=',yt:10:6,'; y2=',y2[1]:10:6,'; y4=',y4[1]:10:6,'; e2=',e2:12:9,'; e4=',e4:12:9,'.');
i_screen:=0;
end;
if t+h>tf-0.000001 then
begin
h:=tf-t;
t:=tf-h;
i_screen:=h_screen;
end;
RK2(t,h,y2);
RK4(t,h,y4);
t:=t+h;
i_screen:=i_screen+h;
end;
writeln;
writeln(' h=',(h):8:5,'; e2max=',e2max:16:8,'; e4max=',e4max:16:8,' ',n:8);
writeln(outfile);
writeln(outfile,' h=',h:6:5,'; e2max=',e2max:10:8,'; e4max=',e4max:10:8,'.');
close(outfile);
readkey;
end.
Функция вычисления точного решения
function clearsolve(t:real):real;
begin
clearsolve:=exp(-t)*(cos(t)+sin(t)+t*sin(t));
end;
Процедура вычисления правых частей системы уравнений в нормальной форме Коши
procedure right(t:real;var x,f:vector_n);
begin
f[1]:=x[2];
f[2]:=2*exp(-t)*cos(t)-2*x[1]-2*x[2];
end;
Процедура RK2
procedure RK2(t:real;h:real;var x:vector_n);
var h4, H33: real;
f1, f2, xr: vector_n;
begin
h4:=0.25*h;
H33:=0.66666667*h;
right(t,x,f1);
for j:=1 to 2 do xr[j]:=x[j]+H33*f1[j];
right(t+H33,xr,f2);
for j:=1 to 2 do x[j]:=x[j]+h4*(f1[j]+3*f2[j]);
end;
Процедура RK4
procedure RK4(t:real;h:real;var x:vector_n);
var H3,h6:real;
f,fs,xr:vector_n;
begin
H3:=0.5*h;
h6:=0.166666667*h;
right(t,x,fs);
for j:=1 to 2 do xr[j]:=x[j]+H3*fs[j];
right(t+H3,xr,f);
for j:=1 to 2 do
begin
xr[j]:=x[j]+H3*f[j];
fs[j]:=fs[j]+2*f[j];
end;
right(t+H3,xr,f);
for j:=1 to 2 do
begin
xr[j]:=x[j]+h*f[j];
fs[j]:=fs[j]+2*f[j];
end;
right(t+h,xr,f);
for j:=1 to 2 do
x[j]:=x[j]+h6*(f[j]+fs[j]);
end;
ПРИЛОЖЕНИЕ Б
(Справочное)
Схема основной программы
Схема функции вычисления точного решения
Схема процедуры вычисления правых частей системы уравнений в нормальной форме Коши
Схема процедуры RK2
Схема процедуры RK4
ПРИЛОЖЕНИЕ В
(Справочное)
Результаты экспериментов
t= 0.000000; yt= 1.000000; y2= 1.000000; y4= 1.000000; e2=0.000000; e4=0.000000
t= 0.400000; yt= 0.982855; y2= 0.983425; y4= 0.982856; e2=0.000570; e4=0.000001
t= 0.800000; yt= 0.893242; y2= 0.893209; y4= 0.893244; e2=0.000033; e4=0.000001
t= 1.200000; yt= 0.726735; y2= 0.725954; y4= 0.726735; e2=0.000781; e4=0.000001
t= 1.600000; yt= 0.518812; y2= 0.517695; y4= 0.518812; e2=0.001117; e4=0.000000
t= 2.000000; yt= 0.312861; y2= 0.311906; y4= 0.312860; e2=0.000955; e4=0.000001
t= 2.400000; yt= 0.141446; y2= 0.140967; y4= 0.141444; e2=0.000479; e4=0.000001
t= 2.800000; yt= 0.020112; y2= 0.020173; y4= 0.020111; e2=0.000062; e4=0.000001
t= 3.200000; yt= 0.050686; y2=-0.050207; y4=-0.050687; e2=0.000480; e4=0.000001
t= 3.600000; yt=-0.080123; y2=-0.079431; y4=-0.080123; e2=0.000692; e4=0.000001
t= 4.000000; yt=-0.081279; y2=-0.080573; y4=-0.081279; e2=0.000705; e4=0.000000
t= 4.400000; yt=-0.066862; y2=-0.066283; y4=-0.066862; e2=0.000579; e4=0.000000
t= 4.800000; yt=-0.046829; y2=-0.046441; y4=-0.046829; e2=0.000388; e4=0.000000
t= 5.200000; yt=-0.027632; y2=-0.027436; y4=-0.027632; e2=0.000196; e4=0.000000
t= 5.600000; yt=-0.012539; y2=-0.012497; y4=-0.012538; e2=0.000042; e4=0.000000
t= 6.000000; yt=-0.002468; y2=-0.002527; y4=-0.002468; e2=0.000058; e4=0.000000
t= 6.400000; yt= 0.003083; y2= 0.002977; y4= 0.003083; e2=0.000106; e4=0.000000
t= 6.800000; yt= 0.005261; y2= 0.005147; y4= 0.005261; e2=0.000113; e4=0.000000
t= 7.200000; yt= 0.005313; y2= 0.005217; y4= 0.005313; e2=0.000096; e4=0.000000
t= 7.600000; yt= 0.004292; y2= 0.004223; y4= 0.004292; e2=0.000068; e4=0.000000
t= 8.000000; yt= 0.002938; y2= 0.002898; y4= 0.002938; e2=0.000040; e4=0.000000
t= 8.400000; yt= 0.001690; y2= 0.001673; y4= 0.001690; e2=0.000017; e4=0.000000
t= 8.800000; yt= 0.000742; y2= 0.000741; y4= 0.000742; e2=0.000001; e4=0.000000
t= 9.200000; yt= 0.000131; y2= 0.000139; y4= 0.000131; e2=0.000008; e4=0.000000
t= 9.600000; yt=-0.000192; y2=-0.000181; y4=-0.000192; e2=0.000011; e4=0.000000
t=10.000000; yt=-0.000310; y2=-0.000299; y4=-0.000310; e2=0.000010; e4=0.000000
h=0.00000; e2max=0.00111851; e4max=0.00000132.