Дифференциальные уравнения
Метод конечных разностей
Конечные разности - самый прямолинейный путь от дифференциального уравнения к коду. Один шаблон Тейлора, одна матрица, один цикл. Именно так считают цены опционов на Wall Street и температуру в ядерных реакторах.
- **Ценообразование опционов:** уравнение Блэка-Шоулза решается схемой Кранка-Николсон на каждом торговом дне. Основа любой системы управления рисками.
- **Термосимуляция:** конечные разности для теплопроводности в GPU/CPU - стандарт проектирования систем охлаждения от Intel до Tesla.
- **Обработка изображений:** лапласиан на сетке - оператор выделения границ. Bilateral filter, guided filter - это дискретизованные PDE.
Предварительные знания
Разностные шаблоны из ряда Тейлора
Производная определена через предел: u'(x) = lim_{h→0} (u(x+h)-u(x))/h. На компьютере предел брать нельзя - но можно остановиться на малом, но конечном h. Какой формулой пользоваться? Ответ даёт ряд Тейлора: разложить u(x±h) в окрестности x, скомбинировать, посмотреть, какие производные сокращаются. Чем больше сокращений - тем выше порядок точности.
Разложение Тейлора: u(x+h) = u(x) + h·u' + h²u''/2 + h³u'''/6 + ... Комбинируя разложения, получаем конечно-разностные аппроксимации: **u'(x) ≈ (u(x+h) - u(x))/h** (вперёд, O(h)); **u'(x) ≈ (u(x+h) - u(x-h))/(2h)** (центральная, O(h²)); **u''(x) ≈ (u(x+h) - 2u(x) + u(x-h))/h²** (вторая, O(h²)).
| Формула | Аппроксимация | Порядок | Шаблон |
|---|---|---|---|
| (u_{i+1}-u_i)/h | u'_i | O(h) | → вперёд |
| (u_{i+1}-u_{i-1})/(2h) | u'_i | O(h²) | ↔ центральная |
| (u_{i+1}-2u_i+u_{i-1})/h² | u''_i | O(h²) | 3-точечная |
| (-u_{i+2}+4u_{i+1}-3u_i)/(2h) | u'_i | O(h²) | → вперёд 2-й |
| (-u_{i+2}+16u_{i+1}-30u_i+16u_{i-1}-u_{i-2})/(12h²) | u''_i | O(h⁴) | 5-точечная |
Центральная разность (u_{i+1}-u_{i-1})/(2h) для u'_i имеет погрешность O(h²). Что это означает при h=0.01?
Схема Кранка-Николсон
Для ∂u/∂t = α∂²u/∂x²: **схема Кранка-Николсон** - среднее явной и неявной схем. Временная производная: (u_i^{n+1} - u_i^n)/Δt. Пространственная: среднее центральных разностей на слоях n и n+1. Результат: **u^{n+1} = u^n + (αΔt/2)(D²u^{n+1} + D²u^n)**, где D² - трёхточечный оператор.
**Схема Кранка-Николсон:** безусловно устойчива для всех Δt (в отличие от явной FTCS: r ≤ 1/2). Точность: O(Δt²) + O(Δx²) - второй порядок по обеим переменным. Цена: на каждом шаге - решение трёхдиагональной СЛАУ.
Кранк, Николсон и Манхэттенский проект (1947)
В 1947 году Джон Кранк и Филлис Николсон (один из первых женщин-вычислителей в Манчестере) опубликовали неявную схему для уравнения теплопроводности. Контекст: Манхэттенский проект и послевоенные расчёты диффузии нейтронов в реакторах. До них использовали явную схему - и при попытке посчитать тепловую динамику топливных стержней с малым Δx схема взрывалась. Кранк и Николсон смешали явный и неявный подходы 50/50 и получили безусловную устойчивость + второй порядок точности. Это позволило считать модели реакторов на доступных компьютерах того времени (Ferranti Mark I в 1951 имел 256 байт памяти).
Сегодня Кранк-Николсон - стандарт для уравнения Блэка-Шоулза на Wall Street и для термических симуляций в EDA-инструментах проектирования чипов.
Почему схема Кранка-Николсон безусловно устойчива, а явная FTCS - условно?
Анализ устойчивости фон Неймана
**Анализ фон Неймана:** предположить, что ошибка разложена в ряд Фурье: ε = Σ Gᵐ·e^{ikx}. Подставить в разностную схему, найти **коэффициент усиления G(k)**. Схема устойчива ⟺ |G(k)| ≤ 1 для всех k. Для FTCS: G = 1 - 4r sin²(kh/2). |G| ≤ 1 ⟺ r ≤ 1/2.
| Схема | G(k) | Условие устойчивости |
|---|---|---|
| Явный Эйлер (FTCS) | 1 - 2r(1-cos kh) | r = αΔt/Δx² ≤ 1/2 |
| Кранк-Николсон | (1-r(1-cos kh))/(1+r(1-cos kh)) | Всегда ≤ 1 |
| Неявный Эйлер | 1/(1+2r sin²(kh/2)) | Всегда ≤ 1 |
| Лакс (волновое) | (cos kh) ± i·c·Δt/Δx·sin kh | CFL: c·Δt/Δx ≤ 1 |
Что означает |G(k)| = 1.001 для некоторого k в анализе фон Неймана?
Приложения: опционы и GPU симуляция
**Уравнение Блэка-Шоулза** для цены европейского опциона: ∂V/∂t + (1/2)σ²S²∂²V/∂S² + rS∂V/∂S - rV = 0. Замена S = e^x переводит к ОДУ с постоянными коэффициентами, которое решается схемой Кранка-Николсон или методом конечных разностей.
**GPU-параллелизм:** разностная схема для PDE - идеальна для GPU. Каждая точка сетки независима (явные схемы) или требует трёхдиагональной солвер (неявные). CuPy = NumPy на GPU. CUDA Thrust - трёхдиагональный Thomas на GPU.
Почему уравнение Блэка-Шоулза решают методом КД, а не Монте-Карло?
Ключевые идеи
- **Разностные шаблоны:** вперёд O(h), центральная O(h²), вторая производная O(h²). Выведены из ряда Тейлора.
- **Кранк-Николсон:** безусловно устойчив, O(Δt²+Δx²). Трёхдиагональная система на каждом шаге.
- **Анализ фон Неймана:** |G(k)| ≤ 1 для всех k - схема устойчива. FTCS: r ≤ 1/2.
- **Блэк-Шоулз:** PDE для цены опциона → КД-схема Кранка-Николсон.
Связанные темы
Метод конечных разностей связывает численные методы с PDE приложениями:
- Метод Эйлера и Рунге-Кутта — FTCS = явный Эйлер по времени + КД по пространству
- Метод конечных элементов — FEM обобщает КД на неструктурированные сетки и слабую формулировку
Вопросы для размышления
- Анализ фон Неймана предполагает постоянные коэффициенты и периодические граничные условия. Для уравнения с переменными коэффициентами этот анализ локальный. Что это означает для практики?
- Уравнение Блэка-Шоулза параболическое. Почему его решают в обратном времени (от T до 0)?
- Конечные разности на GPU: явные схемы параллельны идеально. Почему неявные (трёхдиагональные системы) сложнее распараллелить?