Давление в центре планеты
Сейсмологи по тому, как сквозь планету проходят волны землетрясений, восстанавливают, как её плотность меняется с глубиной, а дальше остальное считают интегралами. Планета состоит из k сферических слоёв, и в каждом слое плотность — многочлен от относительного радиуса x = r/R:
$$\rho(r) = c_0 + c_1 x + c_2 x^2 + c_3 x^3 \ \text{г/см}^3 .$$На границах слоёв плотность обычно прыгает. Найдите три числа:
- массу M = ∫0R 4πr²ρ(r) dr;
- момент инерции в долях MR²: I/(MR²), где I = ∫0R (8π/3)r⁴ρ(r) dr;
- давление в центре p = ∫0R ρ(r) g(r) dr, где g(r) = G·M(r)/r² — сила тяжести на радиусе r, а M(r) — масса внутри этого радиуса; G = 6,6743·10⁻¹¹ м³/(кг·с²).
Третий тест — Земля по модели PREM (Дзевоньский и Андерсон, 1981).
Вход
В первой строке k (1 ≤ k ≤ 100). Дальше k строк от центра наружу: rот rдо c₀ c₁ c₂ c₃ — границы слоя в километрах и коэффициенты плотности. Первый слой начинается в нуле, каждый следующий — там, где кончился предыдущий, последний кончается на радиусе планеты R ≤ 100 000 км. Плотность внутри каждого слоя положительна.
Выход
Три числа через пробел: M в килограммах, I/(MR²) и p в гигапаскалях. Каждое принимается с относительной ошибкой до 10⁻⁶ (для чисел меньше единицы — с абсолютной).
Примеры
Пример 1
1 0 6371 5.5 0 0 0
5.957638043e+24 0.400000000 171.634633287
Пример 2
2 0 3480 11 0 0 0 3480 6371 4.5 0 0 0
6.021896787e+24 0.346521455 330.390819003
Ограничения
- Время · Python
- 5000 мс
- Время · Go
- 1000 мс
- Память
- 256 МБ
- Проверка
- 9 тестов
Разбор
Открывайте по ступеням: каждая следующая подсказывает больше. Можно сразу к решению — но задача, решённая с одной подсказкой, запоминается лучше.
1Подсказка 1
Масса и момент инерции — интегралы многочленов, если брать их по слоям. Что мешает взять их по всему радиусу сразу?
2Подсказка 2
В интеграле давления под знаком стоит M(r) — сам интеграл. Внутри одного слоя он тоже многочлен. Насколько гладкая функция ρ(r)g(r) внутри слоя и что с ней на границе?
3Идея
Все трудности — на границах слоёв: там плотность прыгает, и любая формула с узлами, разбросанными поверх границ, теряет порядок и ошибается пропорционально шагу. Внутри слоя всё гладко. Значит, интегрировать нужно по каждому слою отдельно, с узлами на его концах. Массу и момент инерции слоя дают первообразные: 4πr²ρ(r) — многочлен степени до 5, и ∫ rm = rm+1/(m + 1). Давление — численно, но тоже послойно.
4Разбор
Сначала массы внутри каждой границы: M₀ = 0, Mi+1 = Mi + масса слоя i по первообразной. Последняя — масса планеты. Момент инерции — та же сумма с r⁴ вместо r².
Внутри слоя i масса под радиусом r равна Mi плюс интеграл от нижней границы слоя до r — снова по первообразной, так что g(r) вычисляется точно в любой точке. Подынтегральная функция давления ρ·g внутри слоя гладкая, хотя и не многочлен (в ней деление на r²). Эталон берёт 64 куска на слой и по три точки Гаусса на кусок — ошибка далеко за пределами шестого знака. Подошёл бы и Симпсон с несколькими сотнями кусков на слой: важно не какая формула, а то, что ни один кусок не перекрывает границу.
Проверить себя помогает первый тест: однородная планета. Там I = 0,4MR², а давление в центре 2πGρ²R²/3 — для плотности 5,5 г/см³ и радиуса Земли это 171,6 ГПа. У настоящей Земли — 364 ГПа: тяжёлое ядро в центре давит сильнее.
4+Решение
5Сложность
Масса и момент — O(k). Давление — O(k · m) вычислений, где m — число узлов на слой (здесь 192), каждое с вычислением M(r) по первообразной за O(1).
6Типичные ошибки
- Равный шаг через весь радиус. Даже 20 000 узлов Симпсона вслепую ошибаются в массе на 10⁻⁵ на втором тесте: на каждой границе кусок с узлами по обе стороны скачка даёт ошибку порядка шага. На первом, однородном, тесте этого не видно — скачков нет.
- M(r) накопленными трапециями. Ошибка в массе под радиусом переходит в g и в давление.
- Деление на ноль в центре. g(0) = 0, но формула M(r)/r² там не определена. Узлы Гаусса лежат строго внутри кусков и центр не трогают; у Симпсона значение в нуле придётся задать вручную.
- Единицы. Радиусы во входе в километрах, плотность в г/см³, а G — в метрах и килограммах: 1 г/см³ = 1000 кг/м³.