Численные методы решения обыкновенных дифференциальных уравнений и систем: Методические указания к курсовой работе, страница 8

Программа

Программу для решения уравнения предлагается составить самостоятельно. Рекомендации.

1)  Все параметры уравнения задать с помощью оператора ввода

2)  Правую часть уравнения описать в программе как функцию пользователя

3)  Результаты вычислений вывести в файл

Файл результатов включить в отчет и сравнить полученные результаты с результатами расчетов Excel. На графике представить решения методом Эйлера и методом Рунге –Кутта.

2.2. Расчет колебаний физического маятника

(Пример решения дифференциального уравнения второго порядка)

2.2.1. Постановка и математическая модель задачи

Колебания физического маятника описываются уравнением второго порядка

                                   (2.11)

Здесь    y - угол отклонение маятника от положения равновесия, с - параметр, причем  c = l /g , где l - длина маятника, g - ускорение свободного падения, f(x,y) - вынуждающая сила, действующая на маятник в текущий момент времени. Если колебания являются свободными, т.е. на маятник не действуют внешние силы, то f(x,y)=0

Рассмотрим простейшее уравнение колебаний - уравнение свободных колебаний физического маятника

                                             (2.12)

Пусть в начальный момент времени маятник отклонен от начального положения на угол y₀ . Пусть начальная скорость маятника . Т.е. можно считать, что в начальный момент времени маятник отклонен и отпущен без каких-либо толчков. Мы получили задачу Коши для уравнения второго порядка.

2.2.2. Расчетные формулы

Перейдем от  уравнения (2.12) к системе уравнений с помощью замены :

с начальными условиями y(0)=y₀, z(0)=0

Зададим начальные значения   x=ay=y₀ z=y1₀ .

Следующие значения вычисляются по формулам (1.12-1.13) до достижения переменной  x  граничного значения   x_гр= b+h/2

2.2.3. Алгоритм решения и программа

Алгоритм решения системы ОДУ полностью аналогичен алгоритму решения уравнения первого порядка. Поэтому можно сделать копии блок-схемы и программы предыдущей задачи и внести в них необходимые изменения. Для задания правых частей системы опишем функции  f1(x,y,z)=z,   f2(x,y,z)=-c*sin(y) . Результатом вычислений  является таблица значений величин x, y, где у – искомая функция. Рекомендуется  вывести и значения функции , т. к. хотя в данной задаче искомой является  функция y(x), но z=y’ часто также имеет физический смысл и позволяет больше узнать об изучаемом процессе.

2.2.4. Выполнение расчетов

Уравнение  (2.12) решалось на отрезке  [0, 6.6]  с параметром  c=4,9

В результате расчетов по программе на Паскале получена таблица  2.2.

                                    Таблица 2.2

____________________________

    x      y        z

____________________________

  0.00   0.50000   0.00000

  0.20   0.45369  -0.45624

  0.40   0.32285  -0.83256

  0.60   0.13096  -1.05615

  0.80  -0.08598  -1.07845

  1.00  -0.28641  -0.89435

  1.20  -0.43250  -0.54463

  1.40  -0.49775  -0.09874

  1.60  -0.47082   0.36407

  1.80  -0.35635   0.76355

  2.00  -0.17473   1.02439