The pressure at the centre of a planet
From how earthquake waves pass through a planet, seismologists reconstruct how its density changes with depth, and the rest is computed by integrals. The planet consists of k spherical layers, and in each layer the density is a polynomial in the relative radius x = r/R:
$$\rho(r) = c_0 + c_1 x + c_2 x^2 + c_3 x^3 \ \text{g/cm}^3 .$$At layer boundaries the density usually jumps. Find three numbers:
- the mass M = ∫0R 4πr²ρ(r) dr;
- the moment of inertia as a fraction of MR²: I/(MR²), where I = ∫0R (8π/3)r⁴ρ(r) dr;
- the pressure at the centre p = ∫0R ρ(r) g(r) dr, where g(r) = G·M(r)/r² is gravity at radius r, and M(r) the mass inside that radius; G = 6.6743·10⁻¹¹ m³/(kg·s²).
The third test is the Earth by the PREM model (Dziewoński and Anderson, 1981).
Input
The first line holds k (1 ≤ k ≤ 100). Then k lines from the centre outwards: rfrom rto c₀ c₁ c₂ c₃, the boundaries of the layer in kilometres and the density coefficients. The first layer starts at zero, each next one where the previous ended, and the last ends at the planet’s radius R ≤ 100,000 km. The density inside each layer is positive.
Output
Three numbers separated by spaces: M in kilograms, I/(MR²) and p in gigapascals. Each is accepted with a relative error of up to 10⁻⁶ (for numbers below one, an absolute error).
Examples
Example 1
1 0 6371 5.5 0 0 0
5.957638043e+24 0.400000000 171.634633287
Example 2
2 0 3480 11 0 0 0 3480 6371 4.5 0 0 0
6.021896787e+24 0.346521455 330.390819003
Limits
- Time · Python
- 5000 ms
- Time · Go
- 1000 ms
- Memory
- 256 MB
- Judging
- 9 tests
Walkthrough
Open it step by step: each step gives away a little more. You can jump straight to the solution, but a problem solved with one hint stays with you longer.
1Hint 1
The mass and the moment of inertia are integrals of polynomials if taken layer by layer. What gets in the way of taking them over the whole radius at once?
2Hint 2
The pressure integral has M(r) under the integral sign, itself an integral. Inside one layer it is a polynomial too. How smooth is the function ρ(r)g(r) inside a layer, and what happens to it at a boundary?
3The idea
All the difficulties are at layer boundaries: there the density jumps, and any rule with nodes scattered across the boundaries loses its order and errs in proportion to the step. Inside a layer everything is smooth. So the integration must be done for each layer separately, with nodes at its ends. The mass and moment of inertia of a layer are given by antiderivatives: 4πr²ρ(r) is a polynomial of degree up to 5, and ∫ rm = rm+1/(m + 1). The pressure numerically, but also layer by layer.
4Walkthrough
First the masses inside each boundary: M₀ = 0, Mi+1 = Mi + the mass of layer i by its antiderivative. The last is the mass of the planet. The moment of inertia is the same sum with r⁴ instead of r².
Inside layer i the mass under radius r equals Mi plus the integral from the layer’s lower boundary to r, again by the antiderivative, so g(r) is computed exactly at any point. The integrand of the pressure, ρ·g, is smooth inside a layer, though not a polynomial (it involves division by r²). The reference solution takes 64 pieces per layer and three Gauss points per piece, with an error far beyond the sixth digit. Simpson with a few hundred pieces per layer would also do: what matters is not which rule, but that no piece straddles a boundary.
The first test helps to check yourself: a uniform planet. There I = 0.4MR², and the pressure at the centre is 2πGρ²R²/3, which for a density of 5.5 g/cm³ and the Earth’s radius is 171.6 GPa. The real Earth has 364 GPa: the heavy core at the centre presses harder.
4+Solution
5Complexity
Mass and moment are O(k). The pressure is O(k · m) evaluations, where m is the number of nodes per layer (192 here), each with M(r) computed by the antiderivative in O(1).
6Common mistakes
- An equal step across the whole radius. Even 20,000 Simpson nodes placed blindly err in the mass by 10⁻⁵ on the second test: at each boundary a piece with nodes on both sides of the jump gives an error of the order of the step. On the first, uniform, test this is not visible: there are no jumps.
- M(r) by accumulated trapezoids. An error in the mass under the radius carries over into g and into the pressure.
- Division by zero at the centre. g(0) = 0, but the formula M(r)/r² is undefined there. Gauss nodes lie strictly inside the pieces and do not touch the centre; for Simpson the value at zero has to be set by hand.
- Units. The radii in the input are in kilometres, the density in g/cm³, and G in metres and kilograms: 1 g/cm³ = 1000 kg/m³.