3.2: Чисельні методи

Інженерні спеціальності - це більшість студентів у вигляді курсу фізики, для якого розроблена ця книга, тому, швидше за все, ви потрапляєте в цю категорію. Хоча ви, безумовно, визнаєте, що фізика є важливою частиною вашого навчання, якщо ви мали будь-який вплив на те, як інженери дійсно працюють, ви, ймовірно, скептично ставитеся до смаку вирішення проблем, викладеного в більшості наукових курсів. Ви розумієте, що не дуже багато практичних інженерних розрахунків потрапляють у вузьке коло завдань, для яких точне рішення можна розрахувати за допомогою аркуша паперу і гострого олівця. Реальні проблеми, як правило, складні, і, як правило, їх потрібно вирішувати за допомогою числового хрусту на комп'ютері, хоча ми часто можемо отримати розуміння, працюючи з простими наближеннями, які мають алгебраїчні рішення. Мало того, що числове вирішення проблем є більш корисним у реальному житті, це також навчальний; як початківець студент фізики, я дійсно відчував, що я зрозумів рух снаряда після того, як я працював з ним обома способами, використовуючи алгебру, а потім комп'ютерну програму. (Це було ще в ті часи, коли 64 кілобайт пам'яті вважалося багато.)

У цьому розділі ми почнемо з того, як застосувати числові методи до деяких простих задач, для яких ми знаємо відповідь у «замкнутій формі», тобто єдиного алгебраїчного виразу без будь-якого числення або нескінченних сум. Після цього ми вирішимо проблему, яка зробила б вас всесвітньо відомою, якби ви могли зробити це в сімнадцятому столітті, використовуючи папір і перо! Перш ніж продовжити, вам слід прочитати Додаток 1 на сторінці 908, який знайомить вас з мовою програмування Python.

Спочатку давайте вирішимо тривіальну проблему пошуку того, скільки часу потрібно об'єкту, що рухається зі швидкістю v, щоб проїхати пряму відстань. Ця відповідь закритої форми, звичайно, dist/v, але справа в тому, щоб ввести методи, які ми можемо використовувати для вирішення інших проблем такого типу. Основна ідея полягає в тому, щоб розділити відстань вгору на n рівних частин і скласти час, необхідний для проходження всіх частин. Наступна функція Python виконує цю роботу. Зауважте, що вам не слід вводити номери рядків ліворуч, і вам також не потрібно вводити коментарі. Я опустив підказки >>> і... з метою економії місця.

import math
def time1(dist,v,n):
  x=0 # Initialize the position.
  dx = dist/n # Divide dist into n equal parts.
  t=0 # Initialize the time.
  for i in range(n):
    x = x+dx # Change x.
    dt=dx/v # time=distance/speed
    t=t+dt # Keep track of elapsed time.
  return t

Скільки часу потрібно для переміщення 1 метра при постійній швидкості 1 м/с? Якщо ми це зробимо,

>>> print(time1(1.0,1.0,10)) # dist, v, n
0.99999999999999989

Python дає очікувану відповідь, розділивши відстань на десять рівних відрізків 0,1 метра, і склавши десять 0,1-секундних разів, необхідних для проходження кожного з них. Оскільки об'єкт рухається з постійною швидкістю, навіть не має значення, встановлюємо n на 10, 1 або мільйон:

>>> print(time1(1.0,1.0,1)) # dist, v, n
1.0

Тепер давайте зробимо приклад, де відповідь не очевидна для людей, які не знають обчислення: скільки часу потрібно об'єкту, щоб впасти через висоту h, починаючи з відпочинку? З прикладу 8 на сторінці 83 ми знаємо, що точна відповідь, знайдена за допомогою обчислення, є\(\sqrt{2h/g}\). Давайте подивимося, чи можемо ми відтворити цю відповідь чисельно. Основна відмінність цієї програми від попередньої полягає в тому, що тепер швидкість не є постійною, тому нам потрібно оновлювати її, коли ми йдемо далі. Збереження енергії дає\(mgh=(1/2)mv^2+mgy\) за швидкість\(v\) на висоті\(y\), так що\(v=-\sqrt{2g(h-y)}\). (Ми вибираємо негативний корінь, оскільки об'єкт рухається вниз, а наша система координат має позитивну\(y\) вісь, спрямовану вгору.)

import math
def time2(h,n):
  g=9.8 # gravitational field
  y=h # Initialize the height.
  v=0 # Initialize the velocity.
  dy = -h/n # Divide h into n equal parts.
  t=0 # Initialize the time.
  for i in range(n):
    y = y+dy # Change y. (Note dy<0.)
    v = -math.sqrt(2*g*(h-y)) # from cons. of energy
    dt=dy/v # dy and v are <0, so dt is >0
    t=t+dt # Keep track of elapsed time.
  return t

Для h =1,0 м результат закритої форми дорівнює\(\sqrt{2\cdot1.0\ \text{m}/9.8\ \text{m}/\text{s}^2}=0.45\ \text{s}\). З падінням, розділеним лише на 10 однакових інтервалів висоти, числова техніка забезпечує досить паршиве наближення:

>>> print(time2(1.0,10)) # h, n
0.35864270709233342

Але збільшивши n до десяти тисяч, ми отримуємо відповідь, яка настільки близька, наскільки нам потрібно, враховуючи обмежену точність вихідних даних:

>>> print(time2(1.0,10000)) # h, n
0.44846664060793945

Тонкий момент полягає в тому, що ми змінили y в рядку 9, а потім на лінії 10 ми обчислили v, який залежить від y. Оскільки у змінюється тільки на десятитисячних метра з кожним кроком, ви можете подумати, що це не буде мати великої різниці, і ви б майже праві, за винятком однієї невеликої проблеми: якщо ми поміняли рядки 9 і 10, то в перший раз через цикл, ми б V = 0, що б виробляють помилку поділу на нуль, коли ми обчислювали dt! Насправді, що було б найбільш доцільним було б обчислити швидкість на висоті y і швидкість на висоті y+dy (нагадуючи, що dy негативний), усереднити їх разом, і використовувати це значення y для обчислення найкращої оцінки швидкості між ці два пункти. Оскільки прискорення є постійним у цьому прикладі, ця модифікація призводить до створення програми, яка дає точний результат навіть для n =1:

import math
def time3(h,n):
  g=9.8
  y=h
  v=0
  dy = -h/n
  t=0
  for i in range(n):
    y_old = y
    y = y+dy
    v_old = math.sqrt(2*g*(h-y_old))
    v = math.sqrt(2*g*(h-y))
    v_avg = -(v_old+v)/2.
    dt=dy/v_avg
    t=t+dt
  return t

>>> print(time3(1.0,1)) # h, n
0.45175395145262565

a/ Наближення до кривої брахістохрони з використанням полінома третього порядку (суцільна лінія) та полінома сьомого порядку (пунктирне). Останній лише покращує час на чотири мілісекунди.

Тепер ми готові атакувати проблему, яка кинула виклик кращим умам Європи ще в ті часи, коли не було комп'ютерів. У 1696 році математик Йоганн Бернуллі поставив наступний відомий питання. Починаючи з відпочинку, об'єкт безперешкодно ковзає по кривій, що з'єднує точку\((a,b)\) з точкою\((0,0)\). З усіх можливих форм, які може мати така крива, яка доставить об'єкт до місця призначення за мінімально можливий час, і скільки часу це займає? Оптимальна крива називається брахістохроном, від грецького «короткий час». Рішення проблеми брахістохрона ухилялося від самого Бернуллі, а також Лейбніца, який був одним з винахідників обчислення. Англійський фізик Ісаак Ньютон, однак, не спав пізно ввечері після денної роботи над королівським монетним двором, і, за легендою, виготовив алгебраїчний розчин о четвертій ранку. Потім він опублікував його анонімно, але Бернуллі, як кажуть, зауважив, що коли він прочитав його, він миттєво знав зі стилю, що це Ньютон - він міг «сказати лева від позначки його кігтя».

Замість того, щоб намагатися точного алгебраїчного розв'язку, як це зробив Ньютон, ми отримаємо числовий результат для форми кривої та мінімального\(a\) часу, в особливому випадку\(b\) =1.0 m і =1.0 м Інтуїтивно ми хочемо почати з досить крутого падіння, оскільки будь-яку швидкість ми можемо наростити на початку допоможе нам протягом усього іншого руху. З іншого боку, можна переборщити з цією ідеєю: якщо ми опустимося прямо вниз на всю вертикальну відстань, а потім зробимо поворот під прямим кутом, щоб покрити горизонтальну відстань, отриманий час 0,68 с зовсім трохи довше оптимального результату, причина в тому, що шлях надмірно довгий . Існує нескінченно багато можливих кривих, для яких ми могли б обчислити час, але давайте розглянемо поліноми третього порядку,

\[\begin{equation*} y = c_1x+c_2x^2+c_3x^3 , \end{equation*}\]

де ми вимагаємо для\(c_3=(b-c_1a-c_2a^2)/a^3\) того, щоб крива проходила через точку\((a,b)\). Програма Python, наведена нижче, мало чим відрізняється від того, що ми робили раніше. Функція запитує тільки\(c_1\) і\(c_2\), і обчислює\(c_3\) внутрішньо в рядку 4. Оскільки рух двовимірний, ми повинні обчислити відстань між однією точкою і наступною, використовуючи теорему Піфагора, в рядку 16.

import math
def timeb(a,b,c1,c2,n):
  g=9.8
  c3 = (b-c1*a-c2*a**2)/(a**3)
  x=a
  y=b
  dx = -a/n
  t=0
  for i in range(n):
    y_old = y
    x = x+dx
    y = c1*x+c2*x**2+c3*x**3
    dy = y-y_old
    v_old = math.sqrt(2*g*(b-y_old))
    v = math.sqrt(2*g*(b-y))
    v_avg = (v_old+v)/2.
    ds = math.sqrt(dx**2+dy**2) # Pythagorean thm.
    dt=ds/v_avg
    t=t+dt
  return t

Як перше припущення, ми могли б спробувати пряму діагональну лінію\(y=x\), яка відповідає встановленню\(c_1=1\), а всі інші коефіцієнти до нуля. Результат виходить досить тривалий час:

>>> a=1.
>>> b=1.
>>> n=10000
>>> c1=1.
>>> c2=0.
>>> print(timeb(a,b,c1,c2,n))
0.63887656499994161

Нам дійсно потрібна крива, яка дуже крута праворуч, і лещить ліворуч, тож насправді має сенс спробувати\(y=x^3\):

>>> c1=0.
>>> c2=0.
>>> print(timeb(a,b,c1,c2,n))
0.59458339947087069

Це значне поліпшення, і виявляється всього на соту частку секунди від найкоротшого часу! Можна, хоча і не дуже виховно чи цікаво, знайти кращі наближення до кривої брахістохрони, возившись з коефіцієнтами полінома вручну. Реальний сенс цієї дискусії полягав у тому, щоб навести приклад нетривіальної проблеми, яку можна успішно атакувати за допомогою числових методів. Я знайшов перше наближення, показане на малюнку а,

\[\begin{equation*} y = (0.62)x+(-0.93)x^2+(1.31)x^3 \end{equation*}\]

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

Дописувачі та атрибуція

Template:ContribCrowell