使用Python编写模拟程序，生成假设的脉冲星信号并添加引力波扰动，分析脉冲星信号的到达时间变化，拟合模型参数

本文链接：https://blog.csdn.net/grdd231/article/details/138922922

我们将使用Python编写一个模拟程序，生成假设的脉冲星信号并添加引力波扰动，分析脉冲星信号的到达时间变化，并拟合模型参数。这个过程涉及几个主要步骤：生成脉冲星信号、添加引力波扰动、分析信号的到达时间变化、拟合模型参数。

步骤与代码

1. 安装所需的Python库

pip install numpy scipy matplotlib emcee

2. 编写Python代码

以下是完整的代码示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import emcee

# 定义脉冲星信号的参数
np.random.seed(42)
num_pulsars = 10
total_time = 10.0  # 总观测时间，单位为年
time_intervals = np.linspace(0, total_time, 1000)  # 观测时间点

# 生成假设的脉冲星信号
pulsar_frequencies = np.random.uniform(1, 10, num_pulsars)  # 每个脉冲星的频率
pulsar_phases = np.random.uniform(0, 2*np.pi, num_pulsars)  # 每个脉冲星的初始相位

# 生成脉冲星的到达时间信号
def generate_pulsar_signal(time_intervals, frequency, phase):
    return np.sin(2 * np.pi * frequency * time_intervals + phase)

signals = np.array([generate_pulsar_signal(time_intervals, f, p) for f, p in zip(pulsar_frequencies, pulsar_phases)])

# 添加引力波扰动
def add_gravitational_wave(time_intervals, amplitude, frequency):
    return amplitude * np.sin(2 * np.pi * frequency * time_intervals)

gw_amplitude = 0.1
gw_frequency = 1.0 / total_time  # 假设一个周期为总观测时间
gw_signal = add_gravitational_wave(time_intervals, gw_amplitude, gw_frequency)

# 将引力波扰动添加到每个脉冲星信号中
signals_with_gw = signals + gw_signal

# 分析信号的到达时间变化
# 这里我们假设可以直接观测到引力波的影响（实际情况更复杂）
time_delays = signals_with_gw - signals

# 拟合模型参数：我们尝试拟合出引力波的幅度和频率
def model(time_intervals, amplitude, frequency):
    return amplitude * np.sin(2 * np.pi * frequency * time_intervals)

# 使用曲线拟合
params, params_covariance = curve_fit(model, time_intervals, np.mean(time_delays, axis=0))

# 使用 MCMC 进一步优化参数
ndim, nwalkers = 2, 50
pos = [params + 1e-4*np.random.randn(ndim) for i in range(nwalkers)]

def log_probability(theta, x, y):
    amplitude, frequency = theta
    model_y = model(x, amplitude, frequency)
    sigma2 = 0.01 ** 2
    return -0.5 * np.sum((y - model_y) ** 2 / sigma2 + np.log(sigma2))

sampler = emcee.EnsembleSampler(nwalkers, ndim, log_probability, args=(time_intervals, np.mean(time_delays, axis=0)))
sampler.run_mcmc(pos, 5000, progress=True)

samples = sampler.get_chain(discard=100, thin=15, flat=True)

# 提取拟合参数
amplitude_mcmc, frequency_mcmc = np.percentile(samples, 50, axis=0)

# 绘制结果
plt.figure(figsize=(10, 6))

plt.subplot(2, 1, 1)
plt.title("Pulsar Signal with Gravitational Wave Disturbance")
for i in range(num_pulsars):
    plt.plot(time_intervals, signals_with_gw[i], label=f'Pulsar {i+1}')
plt.xlabel("Time (years)")
plt.ylabel("Signal")

plt.subplot(2, 1, 2)
plt.title("Fitted Gravitational Wave")
plt.plot(time_intervals, np.mean(time_delays, axis=0), label='Observed Time Delays', linestyle='dotted')
plt.plot(time_intervals, model(time_intervals, amplitude_mcmc, frequency_mcmc), label='Fitted GW Signal', linestyle='--')
plt.xlabel("Time (years)")
plt.ylabel("Time Delay")
plt.legend()

plt.tight_layout()
plt.show()

print(f"Fitted Amplitude: {amplitude_mcmc}")
print(f"Fitted Frequency: {frequency_mcmc}")