TP5 Числова роздільна здатність диференціальних рівнянь; ІТ-документація для всіх -

Нас цікавить метод Ейлера для приблизної роздільної здатності важкого гармонічного або затухаючого маятника, а також інші методи.

роздільна

Кілька фазових портретів, що демонструють три можливі режими (субкритичний, критичний, надкритичний):

І приємна анімація ідеального важкого маятника в докритичному режимі:

Документація¶

TP5 у Lycée Lakanal (тема написана Арно Бассоном), на важкому маятнику.

  • Дата: середа 01-13 та четвер 01-14 (2016),
  • Автор: Ліліан Бессон, за курс інформатики для всіх курсів підготовки до МП (http://perso.crans.org/besson/infoMP/),
  • Ліцензія: Ліцензія MIT (http://lbesson.mit-license.org).
TP5. l = 0,3¶

Глобальна змінна \ (l = 30 \) см (довжина стрижня).

Глобальна змінна \ (g = 9,80665 \; м/с ^ 2 \) (прискорення сили тяжіння на поверхні Землі).

Оператор \ (f: Y \ mapsto f (Y) \) (1), записаний як функція Python.

  • Аргументи:
    • Y є одновимірним (двоелементним) масивом numpy і повертає масив numpy однакового розміру. \ (Y = [\ theta, \ dot \ theta] \) так \ (Y_0 = \ theta, Y_1 = \ точка \ theta \),
    • (необов’язково) g - прискорення сили тяжіння на поверхні Землі (\ (g = 9,80665 \; м/с ^ 2 \)),
    • (необов’язково) l - довжина стрижня (\ (l = 30 \) см в додатках).
  • Приклад:

Числово розв’яжіть диференціальне рівняння порядку 2 простого маятника (2) методом Ейлера для \ (t \ в [0, t _] \) з кроком у часі \ (h> 0 \) .

Для використання схеми Ейлера нам потрібно записати це рівняння (звичайне порядку 2, але не лінійне!) У вигляді \ (Y '(t) = f (t, Y) \), де \ (Y = [\ theta, \ dot \ theta] \):

І тому написана схема оновлення Ейлера:

Але в нашому випадку дискретні часи \ (t_i \) розташовані рівномірно, тому \ (\ forall i \ in \, \; t_ - t_ = h \) (часовий крок), тому діаграма спрощена в (3):

  • Оператор \ (f: (t, Y) \ mapsto f (t, Y) \) є векторним (розміром \ (2 \ по 2 \)), і він задається функцією f () .
  • Початкові умови задані theta0 та thetaPoint0 .
  • Повертає масив (список), що містить значення \ (u_c \), обчислене в рази \ (t_i = i \ times h \) (для \ (0 \ leq i \ leq t_/h) \) .
  • Аргументи:
    • h - крок у часі (\ (h> 0 \)),
    • tmax - загальна тривалість цифрового моделювання (\ (t_> 0 \)),
    • theta0 \ (= \ theta_0 \) - значення \ (\ theta (t = 0) \),
    • а thetaPoint0 \ (= \ крапка \) - значення \ (\ крапка \ тета (t = 0) = \ ліворуч (\ frac \ right) (t = 0) \) .
  • Гіпотеза: Ми вирішимо рівняння для \ (t \ geq 0 \) .
  • Приклад (питання I.2.c):
  • Крива \ ((t, \ theta (t)) \), ми вже бачимо, що це не зовсім гармонійне рішення (спостерігається незначне посилення через чисельні помилки, введені наближенням Тейлора, що виправдовує оновлення діаграми Ейлера метод):
  • Фазовий портрет (\ ((\ theta (t), \ dot \ theta (t)) \)), який не є колом, і тому ми бачимо, що метод Ейлера не є точним для цього моделювання:
  • Складність часу: \ (O (n) \), з \ (n = \ lceil t_/h \ rceil \) .
  • Складність в пам’яті (питання I.2.c): константи, плюс масив 64-бітних числових знаків із числами (numpy.float64, тобто 4 байти), розміром точно \ ((n, 2) \), з \ (n = \ lceil t_/h \ rceil \) (розмір масиву numpy можна отримати за допомогою функції numpy.shape ()).
  • Приклад:

\ (\ omega_0 = \ sqrt \) - це правильна пульсація гармонічної системи (докладніше див. цю сторінку).

Намалюйте криву \ (\ theta (t) \), значення кута \ (\ theta \) як функцію часу \ (t \) .

  • Аргументи:
    • listTime - це масив (список або масив numpy), що містить значення \ (t_i \),
    • theta містить значення (приблизні чи ні) \ (\ theta_i = \ theta (t_i) \) .
TP5. portraitPhase (thetas, method = "методом Ейлера") [джерело] ¶

Намалюйте фазовий портрет, тобто похідну \ (\ dot \ theta \) залежно від значення кута \ (\ theta \) .

  • Аргумент:
    • theta - це масив (список або масив numpy) розміром \ ((n, 2) \), який містить значення (приблизні чи ні) \ (Y = [\ theta, \ dot \ theta] \ ) у миттєві випадки \ (t_i \) .
TP5. середня точка (h, tmax, theta0, thetaPoint0) [джерело] ¶

Числово розв’яжіть диференціальне рівняння порядку 2 простого маятника (4) методом середньої точки для \ (t \ in [0, t _] \) з кроком у часі \ (h \) .

Для використання цієї іншої схеми нам також потрібно написати це диференціальне рівняння (звичайне порядку 2, але не лінійне!) У вигляді \ (Y '(t) = f (t, Y) \), де \ (Y = [\ theta, \ dot \ theta] \):

Це той самий оператор \ (f \) (обчислений за допомогою функції f () вище). І отже, схема для оновлення методу середньої точки написана:

Але в нашому випадку дискретні часи \ (t_i \) розташовані рівномірно, тому \ (\ forall i \ in \, \; t_ - t_ = h \) (часовий крок), тому діаграма спрощена в (5):

  • Оператор \ (f: (t, Y) \ mapsto f (t, Y) \) є векторним (розміром \ (2 \ по 2 \)), і він такий самий, як і метод Ейлера (заданий функцією f ()).
  • Початкові умови задані theta0 та thetaPoint0 .
  • Повертає масив (список), що містить значення \ (u_c \), обчислене в рази \ (t_i = i \ times h \) (для \ (0 \ leq i \ leq t_/h) \) .
  • Аргументи:
    • h - крок у часі (\ (h> 0 \)),
    • tmax - загальна тривалість цифрового моделювання (\ (t_> 0 \)),
    • theta0 \ (= \ theta_0 \) - значення \ (\ theta (t = 0) \),
    • а thetaPoint0 \ (= \ крапка \) - значення \ (\ крапка \ тета (t = 0) = \ ліворуч (\ frac \ right) (t = 0) \) .
  • Гіпотеза: Ми вирішимо рівняння для \ (t \ geq 0 \) .
  • Результат: Це порядку 2, отже, більш точний, ніж метод Ейлера, який є лише порядку 1.
  • Приклад (питання I.3.a):
  • Крива \ ((t, \ theta (t)) \), на цей раз вона набагато більше схожа на ідеальне гармонійне рішення (ні затухаюче, ні посилене):
  • Фазовий портрет (\ ((\ theta (t), \ dot \ theta (t)) \)), на якому ми із задоволенням зазначаємо, що метод середньої точки є абсолютно точним для роздільної здатності цього диференціального рівняння (4):
  • Складності: \ (O (n) \) у часі, \ (O (n) \) в пам'яті, з \ (n = \ lceil t_/h \ rceil \) .

Інші посилання на метод середньої точки:

  • Ці слайди (слайд 13),
  • Цей предмет ТП (§ 3),
  • Ці інші слайди,
  • (і нічого у Вікіпедії, але. не соромтеся брати участь у Вікіпедії та створювати чи редагувати сторінку методом midpoint).

Немає часу, \ (h = t_/n> 0 \) (у секундах). Має бути „відносно малим” порівняно з \ (t_ \) .

Оператор \ (g: (Y, t) \ mapsto g (Y, t) \), тут просто дорівнює \ (g (Y, t) = f (Y) \) (вказано в f ()).

Використовуючи scipy.integrate.odeint (див. Його документ, з модуля scipy.integrate), ми отримуємо таку криву та фазовий портрет:

TP5. theta0 = 0¶

Початковий кут \ (\ theta_0 = \ theta (0) \) (у радіанах).

TP5. thetaPoint0 = 6¶

Початкова кутова швидкість \ (\ dot = \ dot \ theta (0) \) (у радіанах за секунду).

Максимальний час \ (t_ \) (у секундах).

TP5. nbPoints = 2000¶

Кількість підрахованих балів \ (n \). Для правильного графіку має бути близько тисячі.

TP5. method = "за методом Ейлера" ¶

Вибір методу, щоб змінити заголовок графіки.

Використовуйте scipy.integrate.odeint (), щоб відобразити фазовий портрет простого маятника з декількома траєкторіями.

  • Для цього візьмемо \ (t_ = 2 с \), \ (\ theta_0 = 0 \) та \ (\ dot = k \; \ text/s \), за \ (k \ у [1. 15 ] \) .

Примітка щодо трьох режимів (субкритичний, надкритичний, критичний).

  • Коли \ (\ theta_0 = 0 \) досить малий (\ (\ точка), маятник коливається навколо стабільного положення рівноваги \ (\ theta_ = 0 \) (субкритичний режим),
  • Коли \ (\ theta_0 = 0 \) досить великий (\ (\ dot> \ dot _> \)), він обертається навколо точки прикріплення (надкритичний режим),
  • В критичний режим (що відбувається, коли \ (\ theta_0 = 0 = 2 \ omega_0 \)), маятник нескінченно повільно сходиться до нестабільного положення рівноваги (\ (\ theta = \ pi \)). Цей останній випадок є теоретичним, він ніколи не спостерігається на практиці.

  • У наших прикладах \ (2 \ omega_0 \ simeq 11.4 \), тому, якщо ми виберемо цілі значення від 8 до 15 (включно), 4 буде знаходитися в докритичному режимі (8,9,10,11), тобто 4 кола або еліпси, а 4 будуть у надкритичному режимі (12,13,14,15), тобто 4 синусоїдальних траєкторії на фазовому портреті.
  • На додаток до цього графіка, я додав критичну швидкість, щоб добре її візуалізувати (темно-зеленим). Ми бачимо, що маятник сходиться до (нестійкої) рівноваги, де \ (\ theta = \ pi \):
TP5. період1 (год, тетас) [джерело] ¶

Обчислює приблизне значення періоду \ (T \) динамічної системи, прийнятого в докритичному режимі.

З цією метою ми спостерігатимемо, що \ (\ theta (T/2) = 0 \) та \ (\ theta (t)> 0 \) для \ (t \ in] 0, T/2 [\). Тоді буде достатньо визначити приблизно найменший \ (t> 0 \) такий, що \ (\ theta (t) = 0 \):

Обчислює приблизне значення періоду \ (T \) динамічної системи, прийнятого в докритичному режимі.

Ми обчислюємо \ (\ omega_0 \), потім \ (\ cos (\ theta _) \) за (6), що ми повертаємо назад, щоб мати \ (\ theta_ \) максимальну амплітуду коливань:

Потім обчисліть \ (T \) за цією формулою (7), використовуючи узагальнений інтеграл:

Ми будемо використовувати функцію scipy.integrate.quad () (див. Її документацію?), З того самого модуля scipy.integrate, котрі тут є деякі пояснення). Ця функція quad () є дуже загальним способом обчислення інтегралу чисельно. Для реальних чи багатовимірних функцій, безперервних чи ні, quad () працює дуже добре (для певних інтегралів, але не лише для обмежених областей).

Дивіться цю статтю Вікіпедії про гауссові квадратурні методи, щоб отримати докладнішу інформацію, якщо вам цікаво, або цю сторінку англійською мовою, що виникла в результаті виправлення проекту, який я дав своїм індійським студентам минулого року (квітень 2015).

Ми навіть можемо інтегрувати цими методами в більш складні домени, ніж прості гіпер-прямокутники, див. Ці приклади там (ідея дуже дурна, але дуже розумна!).

Обчислює (і відображає, якщо show = True) анімацію руху важкого маятника.

  • режим може бути "субкритичним", "критичним", "надкритичним" або "нульовим" (\ (\ dot = 0 \), отже, руху немає).
  • show = False за замовчуванням, але анімація відображається перед збереженням, якщо True .

Приклад субкритичний режим (\ (\ точка), з вічними коливаннями та кутом, завжди обмеженим (\ (0 \ leq \ theta \ leq \ theta_) і нижче горизонталі:

  • Приклад майже критична дієта, незважаючи на вибір \ (\ dot = \ dot_> = 2 \ omega_0 \), маятник не досягає нестабільної рівноваги \ (\ theta = \ pi \) (при першому коливанні ми майже в це віримо, але ні):
  • Приклад надкритичний режим (\ (\ точка> \ точка _> \)), з постійними коливаннями, завжди в одному напрямку, \ (0 \ leq \ theta \ до + \ infty \), що розходиться (за межами \ (2 \ pi \), він підраховує повороти):
  • Приклад "нульового" режиму з \ (\ theta_0 = 0 \) та \ (\ dot = 0 \), немає руху (і тому це найбільш захоплююча анімація дня):
TP5. botafumeiro () [джерело] ¶

Тепер ми розглянемо збуджений маятник.

І відео, на якому показано Ботафумейро в дії (через YouTube):