Рівняння Хвилі у Двох Вимірах: Кінцеві Різниці, Рейліха та Стаціонарні Хвилі

Одна рівняння для кожної хвилі

Звук у повітрі, бризки на шкірі барабана, світло вакуумі, струна, що грається, сейсмічна хвиля через скелю — все підкоряється одній другій степені частковій диференціальній рівнянню. У двох вимірах:

∂²u/∂t² = c² ( ∂²u/∂x² + ∂²u/∂y² ) = c² ∇²u ν(x, y, t) = поле зміщення (або тиску) c = швидкість поширення середовища Фізичне читання є прямим. ∇²u — це кривина поля — кількість, наскільки точка відрізняється від середнього свого оточення. Рівняння говорить: точка, яка знаходиться нижче середнього своїх сусідів, прискорюється вгору, а точка, що знаходиться вище середнього, тягнеться вниз. Це закон Гука, поширений на нескінченність, і все хвилеподібне слідує з нього.

∂²u/∂t² = c² ( ∂²u/∂x² + ∂²u/∂y² )  =  c² ∇²u

u(x, y, t) = the displacement (or pressure) field
c          = the propagation speed of the medium

Розбиття на сітку: дві лінії кінцевих різниць

Розмістіть поле на сітці з кроком Δx та інтервалом часу Δt. Замініть обидва другий-порядкового похідні стандартною центральною різницею, яка є точною другого порядку:

∂²u/∂t² ≈ ( u_i,j^(n+1) − 2u_i,j^n + u_i,j^(n−1) ) / Δt² ∇²u ≈ ( u_(i+1),j + u_(i−1),j + u_i,(j+1) + u_i,(j−1) − 4·u_i,j ) / Δx² // п'ятиточковий шаблон Підставте, розв’яжіть для єдиної невідомої — поля на наступному кроці часу — і весь розв’язувач випадає. Це повністю явний метод: без матриць, без лінійної системи, без ітерацій. Зауважте, що це схема третього рівня – оскільки рівняння є другого порядку в часі: потрібно зберігати два попередніх кадри.

const C2 = (c * dt / dx) ** 2; // число Рейнольдса, квадратне for (let j = 1; j < h - 1; j++) for (let i = 1; i < w - 1; i++) { const k = j * w + i; const lap = cur[k - 1] + cur[k + 1] + cur[k - w] + cur[k + w] - 4 * cur[k]; next[k] = 2 * cur[k] - prev[k] + C2 * lap; } [prev, cur, next] = [cur, next, prev]; // обертайте три буфери Це повний робочий 2D розв’язувач хвилі в дев'ять рядків. Ця схема часто називається FDTD — кінцева різниця по часу області — така сама родина методів, що використовуються для рівнянь Максвелла в електромагнетизмі, застосована тут до скалярної хвильової рівняння.

∂²u/∂t² ≈ ( u_i,j^(n+1) − 2u_i,j^n + u_i,j^(n−1) ) / Δt²

∇²u     ≈ ( u_(i+1),j + u_(i−1),j + u_i,(j+1) + u_i,(j−1)
            − 4·u_i,j ) / Δx²            // the 5-point stencil

Хвильовий рівняння у двох вимірах

Явне розв’язання є лише обставинно стабільним. Інформація в чисельній схемі рухається точно на одну клітинку за часний крок (шаблон досягає лише сусідніх клітинок); фізична інформація рухається швидкістю c·Δt. Якщо фізичний хід переганяє шаблон, схема не може його представити, і помилка не просто зростає — вона росте без обмежень. Умова CFL (Courant–Friedrichs–Lewy, 1928) визначає цю вимогу точно. З числом Courant C = c·Δt/Δx, аналіз стабільності von Neumann цієї схеми в d вимірах дає:

C = c·Δt/Δx ≤ 1 / sqrt(d) 1D C ≤ 1 (і C = 1 є точним — "магічний часний крок") 2D C ≤ 1/√2 ≈ 0.7071 3D C ≤ 1/√3 ≈ 0.5774 Порушення цього призводить до того, що симуляція не деградує граціозно: через кілька десятків кроків поле насичується перевіркою ±∞ і кожне значення стає NaN. Якщо ваша хвильова симуляція вибухає, число Courant — це перше, що потрібно обчислити — воно майже завжди є причиною. Корисною межею безпеки є запуск при C ≈ 0.5 в 2D замість теоретичного верхнього значення.

1D дивоцінь вартий знання: точно при C = 1 дискретна схема відтворює точне рішення безперервного рівняння, оскільки швидкість чисельного та фізичного розповсюдження збігається точно. Це єдичний випадок, коли кінцено-різницева схема є точною, і це також попередження — воно не узагальнюється до 2D, де C = 1 комфортно нестабільний.

C = c·Δt/Δx ≤ 1 / sqrt(d)

1D   C ≤ 1        (and C = 1 is exact — the "magic time step")
2D   C ≤ 1/√2  ≈ 0.7071
3D   C ≤ 1/√3  ≈ 0.5774

Численна дисперсія: хвиля, яка прибуває не вчасно

Навіть коли схема стабільна, вона не є ідеальною. Підставляючи плоску хвилю у дискретні рівняння, отримуємо чисельну дисперсійну відношення, яке не зовсім відповідає фізичній: на сітки швидкість розповсюдження залежить від довжини хвилі та напрямку руху. Короткі довжини хвиль – ті, що вирішені лише кількома сітковими точками, рухаються помітно повільніше, ніж довгі, і рухаються трохи різними швидкостями вздовж осей, ніж по діагоналях.

Найбільш помітна ознака – це різкий імпульс, який розвивається з затримкою, розповсюджуючись, навіть якщо справжнє рівняння не є дисперсійним, і колоподібний фронт хвилі, що стає дещо квадратним. Лікарство – це роздільна здатність: забезпечуйте щонайменше 10–20 сіткових точок на довжину хвилі найкоротшої хвилі, яку ви хочете врахувати, і помилка стає незначною. Вищі порядки стінок (9-точковий Лапласіан) зменшують анізотропію для тієї ж сітки

Границі: відбиття, свобода чи зникнення

Що ви вводите в крайні комірки визначає фізику межі, і є три корисних варіанти:

Дірхлетова u = 0 на краю. ЗАМКНЕНІЙ край — барабанна шкіра, прикріплена до її краю, струна, зав'язана. Хвилі відбиваються з інверсією (пік повертається як низь). Нейманів ∂u/∂n = 0 на краю, реалізований шляхом копіювання прилежної внутрішньої комірки в привид комірку. СВОБОДНА межа — відкритий орган. Хвилі відбиваються БЕЗ інверсії. Поглинання Край повинен поводитися ніби домен триває (Мура, PML) назавжди. Застосовується перший порядок Мурової умови 1D одностороннього хвильового рівняння на краю; Ідеально Відповідний шар оточує домен штучно втратною областю, яка майже не відбиває хвилі, і є стандартним у серйозних кодах. Не робити нічого взагалі — залишити крайні комірки зі своїм початковим значенням — це Дірхлетова умова випадково, тому що незакінчений розв'язувач хвиль завжди показує коробку з відбитками. І саме відбиття створює стоячі хвилі та нормальні режими обмеженого домену: на прямокутнику дискретний набір частот точно підходить до межі, і поле стабілізується в стаціонарних шаблонах із фіксованими нодовими лініями. Посипте пісок на вібруючий пластинчастий, і він накопичується вздовж цих ліній — фігури Хладні — який є тим самим явищем, яке відтворює цей розв'язувач на сітці.

Dirichlet   u = 0 at the edge.  A CLAMPED boundary — a drum skin
            nailed to its rim, a string tied down. Waves reflect
            with an inversion (a crest returns as a trough).

Neumann     ∂u/∂n = 0 at the edge, implemented by copying the
            adjacent interior cell into a ghost cell. A FREE
            boundary — an open organ pipe. Waves reflect WITHOUT
            inversion.

Absorbing   The edge should behave as if the domain continued
(Mur, PML)  forever. A first-order Mur condition applies the 1D
            one-way wave equation at the edge; a Perfectly Matched
            Layer surrounds the domain with an artificially lossy
            region that swallows outgoing waves with almost no
            reflection, and is the standard in serious codes.

Часті запитання

Спробуйте наживо

Усе, що вище, працює прямо у вашому браузері — відкрийте Wave Equation і змінюйте параметри під час роботи. Нічого не встановлюється, нічого не завантажується на сервер, уся модель живе в одній вкладці.