用mathematica模擬引力波還是比較複雜的，因為引力波它涉及到兩個黑洞的碰撞併合，是多體問題，需要做很多近似才可以求解，它的模擬計算程式一般用c語言或者c#來寫。具體的話，你可以參考一下北京師範大學曹周鍵老師等人的研究工作。

而且你需要理解EOBNR理論模型，請參考如下工作：

EOB formalism: improved PN [Buonanno and Damour, PRD 1999]

EOBNR: set more parameters and fit NR [Buonanno, 潘奕，等, PRD 2007]

calibrate to NR results [Buonanno, 潘奕，等, PRD 2009 (nonspinning, equal mass black holes), 2010 (nonprecessing, spinning, equal-mass black holes), 2014 (generic mass ratios and aligned spins), 2014 (spinning, precessing)]

不過我建議你從基礎的mathematica程式做起，比如先學會用mathematica計算愛因斯坦引力張量與黎曼張量。大致程式如下，不過我只能給你提供大致的程式程式碼（來自黃宇傲天，他與我一起寫過相對論的科普書）：

先寫出度量張量，比如史瓦西度量

g := {{-(1 - (2 M)/r), 0, 0, 0}, {0, (1 - (2 M)/r)^-1, 0, 0}, {0, 0,

r^2, 0}, {0, 0, 0, r^2 Sin[\[Theta]]^2}};

然後計算聯絡一形式：

Af[i_] := 1/2*IvG.(Diff[g, i] - Pt[i]);

Subscript[\[CapitalGamma]^\[Mu], i\[Nu]]*)

然後計算計算曲率2形式矩陣：

Wf[i_, j_] := Diff[Af[j], i] - Diff[Af[i], j] + Af[i].Af[j] - Af[j].Af[ i];

接下來就是計算各個曲率張量

Ricci[i_, j_] := Sum[Wf[k, j][[k, i]], {k, 1, dim}];(*計算裡奇張量*)

RicMat := Table[Ricci[i, j], {i, 1, dim}, {j, 1, dim}];(*將裡奇張量以矩陣形式儲存*)

RcS := Sum[

IvG[[i, j]]*Ricci[i, j], {i, 1, dim}, {j, 1, dim}];(*計算曲率標量*)

Ein := RicMat - 1/2*RcS*g(*計算愛因斯坦張量，以矩陣形式儲存*)

RicMat := Table[Ricci[i, j], {i, 1, dim}, {j, 1, dim}](*將裡奇張量以矩陣形式儲存*);

這個東西沒有兩年是做不好的，希望以上資訊能幫助到你。如果不是專業人士，基本可以忽略這個回答了，因為很多人估計不理解這裡面的基本邏輯是什麼。不過看你提這樣的問題，我覺得你也許是瞭解相對論吧。

關於LIGO 這裡有一篇報道（https://www.prnewswire.com/news-releases/scientists-pioneer-use-of-deep-learning-for-real-time-gravitational-wave-discovery-300588346.html） Mathematica 構建了一個深度學習框架和 Einstein Toolkit（https://einsteintoolkit.org/） 在超算上面。工程量和資料量都不小。如果只做LIGO 資料分析的話，SE上面有個討論不錯（https://mathematica.stackexchange.com/questions/107874/extracting-and-displaying-data-from-gw150914-hdf5-files）