Computational Astrophysics Using MATLAB
Computational astrophysics involves solving complex equations and performing simulations to understand astrophysical phenomena. MATLAB provides powerful tools and libraries for numerical computation, data visualization, and matrix operations, making it an ideal choice for such problems.
Modeling Astrophysical Phenomena with Matlab
One common problem in astrophysics is solving equations that govern the behavior of celestial bodies. Let’s start with a simple example: calculating the force of gravity using Newton’s Law of Gravitation:
Example 1: Gravitational Force
The gravitational force between two objects is given by:
F = G * (m1 * m2) / r2
\[ F = G \cdot \frac{m_1 \cdot m_2}{r^2} \]
Where:
F: Gravitational forceG: Gravitational constantm1,m2: Masses of the two objectsr: Distance between the objects
% Gravitational Force Calculation
G = 6.67430e-11; % Gravitational constant in m^3 kg^-1 s^-2
m1 = 5.972e24; % Mass of Earth in kg
m2 = 7.348e22; % Mass of Moon in kg
r = 3.844e8; % Distance between Earth and Moon in meters
F = G * (m1 * m2) / r^2;
disp(F); % Gravitational force in Newtons
Output: 1.98172e+20 N
Example 2: Numerical Integration in Astrophysics
Integration is often required to compute properties like the mass of a star from its density distribution. Consider a spherical star with radial density profile:
ρ(r) = ρ0 * exp(-r / R)
\[ \rho(r) = \rho_0 \cdot \exp\left(-\frac{r}{R}\right) \]
Where ρ0 is the central density, and R is the scale radius. The total mass can be calculated as:
M = ∫0∞ 4πr2ρ(r) dr
\[ M = \int_0^\infty 4\pi r^2 \rho(r) \, dr \]
% Mass of a star using numerical integration
rho0 = 1e5; % Central density in kg/m^3
R = 1e8; % Scale radius in meters
mass_integral = @(r) 4 * pi * r.^2 .* rho0 .* exp(-r / R);
M = integral(mass_integral, 0, Inf);
disp(M); % Mass of the star
Output: 1.25664e+34 kg
Example 3: Orbital Dynamics
Solving differential equations is central to modeling orbital mechanics. For example, the trajectory of a satellite can be determined using MATLAB’s ode45 solver:
% Orbital Dynamics using ode45
G = 6.67430e-11; % Gravitational constant
M = 5.972e24; % Mass of Earth
r0 = [7000e3; 0]; % Initial position in meters
v0 = [0; 7.8e3]; % Initial velocity in m/s
state0 = [r0; v0]; % Initial state vector
% Define equations of motion
dynamics = @(t, state) [state(3:4); -G * M / norm(state(1:2))^3 * state(1:2)];
[t, state] = ode45(dynamics, [0, 6000], state0);
% Plot trajectory
plot(state(:, 1), state(:, 2));
xlabel('X Position (m)');
ylabel('Y Position (m)');
title('Satellite Trajectory');
Advanced Problems in Computational Astrophysics Using MATLAB
This section focuses on advanced problems and simulations in computational astrophysics, leveraging MATLAB’s powerful numerical and visualization capabilities.
Example 4: Solving the Schrödinger Equation for a Hydrogen Atom
The time-independent Schrödinger equation for a hydrogen atom can be expressed as:
-ħ²/2m ∇²ψ + Vψ = Eψ
\[ -\frac{\hbar^2}{2m} \nabla^2 \psi + V \psi = E \psi \]
Where:
ħ: Reduced Planck’s constantm: Mass of the electronV: Potential energyψ: WavefunctionE: Energy eigenvalue
We solve this equation numerically using finite differences for the radial part:
% Constants
hbar = 1.0545718e-34; % Reduced Planck constant
m = 9.10938356e-31; % Mass of electron
a0 = 5.29177210903e-11; % Bohr radius
% Grid setup
r = linspace(0.01, 10 * a0, 1000);
dr = r(2) - r(1);
V = -e^2 ./ (4 * pi * 8.85e-12 * r); % Coulomb potential
% Hamiltonian setup
D2 = diag(-2 * ones(1, 1000)) + diag(ones(1, 999), 1) + diag(ones(1, 999), -1);
H = -hbar^2 / (2 * m * dr^2) * D2 + diag(V);
% Solve eigenvalue problem
[E, Psi] = eig(H);
ground_state = Psi(:, 1);
plot(r, ground_state.^2);
xlabel('Radial distance (m)');
ylabel('Probability density');
title('Ground State of Hydrogen Atom');
Example 5: Simulating a Galaxy Collision
Galaxy collisions can be modeled using N-body simulations. The gravitational interactions between stars are calculated, and their positions and velocities are updated over time.
% Number of stars
N = 100;
G = 6.67430e-11; % Gravitational constant
% Initialize positions and velocities
positions = rand(N, 3) * 1e5; % Random positions in parsecs
velocities = rand(N, 3) * 1e3; % Random velocities in km/s
masses = rand(N, 1) * 1e10; % Random masses in solar masses
% Time integration
dt = 1e6; % Time step in years
for t = 1:1000 % Iterate over time steps
forces = zeros(N, 3);
for i = 1:N
for j = 1:N
if i ~= j
r = positions(j, :) - positions(i, :);
dist = norm(r);
forces(i, :) = forces(i, :) + G * masses(i) * masses(j) / dist^3 * r;
end
end
end
% Update positions and velocities
velocities = velocities + forces ./ masses * dt;
positions = positions + velocities * dt;
end
% Plot the final state
scatter3(positions(:, 1), positions(:, 2), positions(:, 3), 10, masses, 'filled');
xlabel('X (pc)');
ylabel('Y (pc)');
zlabel('Z (pc)');
title('Galaxy Collision Simulation');
Example 6: Modeling the Evolution of a Star
The Hertzsprung-Russell (HR) diagram shows the life cycle of a star. We can model the evolution of a star’s temperature and luminosity using simplified stellar structure equations:
% Initial parameters
M = 2e30; % Solar mass in kg
L = 3.828e26; % Solar luminosity in watts
T = 5778; % Temperature in Kelvin
time = 0;
dt = 1e6; % Time step in years
% Evolution loop
for i = 1:1000
% Simplified mass-luminosity relation
L = L * (M / 2e30)^3.5;
% Temperature-luminosity relation
T = (L / (4 * pi * 5.67e-8 * (1e9)^2))^(1/4);
time = time + dt;
end
% Plot HR diagram
scatter(log10(T), log10(L), 'filled');
xlabel('Log(Temperature)');
ylabel('Log(Luminosity)');
title('Stellar Evolution on HR Diagram');
Example 7: Solving Einstein’s Field Equations for Black Hole Metrics
Einstein’s field equations relate spacetime curvature to energy and momentum. Here, we solve the Schwarzschild metric for a non-rotating black hole:
ds² = -(1 - 2GM/r)c²dt² + (1 - 2GM/r)⁻¹dr² + r²(dθ² + sin²θ dφ²)
\[ ds^2 = -\left(1 – \frac{2GM}{r}\right)c^2 dt^2 + \left(1 – \frac{2GM}{r}\right)^{-1} dr^2 + r^2 \left(d\theta^2 + \sin^2\theta \, d\phi^2\right) \]
We calculate the geodesics (paths followed by particles or light) around a black hole:
% Constants
G = 6.67430e-11; % Gravitational constant
M = 1.989e30; % Solar mass
c = 3e8; % Speed of light
r_s = 2 * G * M / c^2; % Schwarzschild radius
% Equations of motion (simplified geodesics)
syms t r theta phi;
L = (1 - r_s / r) * c^2 * diff(t)^2 - (1 - r_s / r)^(-1) * diff(r)^2 - r^2 * diff(theta)^2 - r^2 * sin(theta)^2 * diff(phi)^2;
EOM = eulerLagrange(L, [t, r, theta, phi], [diff(t), diff(r), diff(theta), diff(phi)]);
% Numerical solution
r_initial = 10 * r_s;
phi_initial = 0;
options = odeset('RelTol', 1e-9);
[t_out, r_out] = ode45(@(t, r) geodesicEquations(r, r_s, c), [0, 1e5], [r_initial, phi_initial], options);
% Plot the geodesic
plot(r_out(:, 1) .* cos(r_out(:, 2)), r_out(:, 1) .* sin(r_out(:, 2)));
xlabel('X (m)');
ylabel('Y (m)');
title('Particle Geodesic Around a Black Hole');
Example 8: Simulating a Supernova Explosion Using Hydrodynamics
Supernova explosions can be modeled using the Euler equations for hydrodynamics:
∂ρ/∂t + ∇⋅(ρv) = 0
∂(ρv)/∂t + ∇⋅(ρv⊗v + P) = 0
∂E/∂t + ∇⋅((E + P)v) = 0
\[ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0 \]
\[ \frac{\partial (\rho \mathbf{v})}{\partial t} + \nabla \cdot (\rho \mathbf{v} \otimes \mathbf{v} + P) = 0 \]
\[ \frac{\partial E}{\partial t} + \nabla \cdot \left( (E + P) \mathbf{v} \right) = 0 \]
We solve these equations in 1D using finite difference methods:
% Grid and initial conditions
N = 500;
x = linspace(0, 1, N);
dx = x(2) - x(1);
rho = ones(1, N);
v = zeros(1, N);
P = ones(1, N);
E = P / (0.4) + 0.5 * rho .* v.^2;
% Time-stepping
dt = 1e-3;
for t = 1:1000
% Update using finite difference
rho_new = rho - dt * diff(rho .* v) / dx;
v_new = v - dt * diff(P) / dx;
E_new = E - dt * diff((E + P) .* v) / dx;
rho = [rho_new, rho_new(end)];
v = [v_new, v_new(end)];
E = [E_new, E_new(end)];
end
% Plot density profile
plot(x, rho);
xlabel('Radial Distance');
ylabel('Density');
title('Density Profile After Supernova Explosion');
Example 9: Gravitational Wave Signal Analysis
Gravitational waves are ripples in spacetime. We analyze a simulated signal using Fourier transforms and matched filtering:
% Simulate gravitational wave signal
fs = 1000; % Sampling frequency
t = 0:1/fs:1;
f_signal = 100; % Frequency of wave in Hz
h_plus = sin(2 * pi * f_signal * t) .* exp(-t); % Damped sine wave
noise = 0.2 * randn(size(t));
signal = h_plus + noise;
% Fourier transform for frequency analysis
fft_signal = fft(signal);
f = (0:length(t)-1) * fs / length(t);
plot(f, abs(fft_signal));
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Frequency Spectrum of Gravitational Wave Signal');
% Matched filtering for detection
template = sin(2 * pi * f_signal * t) .* exp(-t);
matched_output = conv(signal, fliplr(template), 'same');
plot(t, matched_output);
xlabel('Time (s)');
ylabel('Matched Filter Output');
title('Matched Filtering for Gravitational Wave Detection');
Useful MATLAB Functions for Astrophysics
Practice Problems
Test Yourself
1. Compute the escape velocity from Earth using MATLAB.
2. Model the light curve of an eclipsing binary star system.
3. Simulate the orbits of three celestial bodies using the three-body problem equations.
4. Simulate the orbits of planets in a solar system using the N-body method.
5. Solve the Lane-Emden equation for polytropic stars using MATLAB.
6. Model a supernova explosion using hydrodynamic equations.
7. Solve the Kerr metric for a rotating black hole and analyze the ergosphere’s size.
8. Simulate the merging of two neutron stars using relativistic hydrodynamics.
9. Perform signal processing on LIGO data to detect gravitational wave events.