该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
Clear["Global`*"]
c = 299792458*10^2(*光速，单位cm/s*)
G = 6.67259*10^-8(*gravitational constant，引力常数，单位cm^3/g*s^2*)
Msun = 1.9891*10^33(*Subscript[M, \[CircleDot]]，太阳质量，单位g*)
Itilder = 0.283(*Overscript[I, ~]，(7)式下面*)
Jtilder = 1.81*10^-2(*Overscript[J, ~]，(23)式下面*)
Mbi = 1.4*Msun(*Subscript[M, b,i]，(29)式下面*)
\[Eta] = 0.01(*\[Eta]，fig1的不同情况*)
t0 = 3000(*Subscript[t, 0]，fig1的不同情况*)
B = 2*10^14(*(29)式下面，单位G*)
R115 = 1(*Subscript[R, 11.5]*)
R = 11.5*10^5(*单位cm，R=11.5km，(29)式下面*)
M14 = 1(*Subscript[M, 1.4]*)
B15 = 1(*Subscript[B, 15]*)
Bt14 = 1(*Subscript[B, t,14]*)
T9 = 1(*Subscript[T, 9]*)
u = B*R^3(*\[Mu]，磁偶极矩，(18)式下面*)
Medot[t_] :=
10^-3*\[Eta]*t^(1/2)*Msun(*(13)式，Subscript[Overscript[M, .], early]*)
Mldot[t_] :=
10^-3*\[Eta]*t0^(13/6)*t^(-5/3)*
Msun(*(14)式，Subscript[Overscript[M, .], late]*)
Mdot[t_] := (1/Medot[t] + 1/Mldot[t])^-1(*(12)式，Overscript[M, .]*)
M[t_] := Mb[t]*(1 + (3*G*Mb[t])/(5*R*c^2))^-1(*(16)式，M*)
Mbdot[t_] :=
Mdot[t]/((1 + 3/5*G*Mb[t]/(R*c^2))^-1 -
Mb[t]*3/5*G/(R*c^2)*(1 + 3/5*G*Mb[t]/(R*c^2))^-2)
II[t_] := Itilder*M[t]*R^2(*I，(6)式下面*)
v3[t_] := \[CapitalOmega][t]/(2*Pi*10^3)(*Subscript[\[Nu], 3]，(2)式下面*)
CC[t_] := 2.3*10^-4*Bt14^2*M14^-1*v3[t]^-2(*C，改写为CC(7)式下面*)
rm[t_] := (u^4/(G*M[t]*Mbdot[t]^2))^(
1/7)(*Subscript[r, m]，(18)式，单位km*)
\[Omega][t_] := \[CapitalOmega][t]/Sqrt[G*M/rm[t]^3](*\[Omega]，(19)式*)
n\[Omega][t_] := 1 - \[Omega][t](*n(\[Omega])，(21)式下面*)
rc[t_] := 16.5*M14^(1/3)*v3[t]^(-2/3)(*Subscript[r, c]，(17)式，单位km*)
tsv = 1.4*10^8*M14*R115^-1*T9^(5/3)(*Subscript[t, sv]，(3)式，单位s*)
tgw[t_] := -24*R115^-4*M14^-1*v3[t]^-6(*Subscript[t, gw]，(2)式，单位s*)
tBt[t_] :=
5.8*R115^-1*M14*B15^-1*Bt14^-1*v3[t](*Subscript[t, B,t]，(5)式，单位s*)
tBgw[t_] :=
3.8*10^14*M14^-1*R115^-2*Bt14^-4*
v3[t]^-4(*Subscript[t, B,gw]，(8)式，单位s*)
Nacc[t_] := n\[Omega][t]*u^2/(rm[t]^3)(*(20)式*)
tdip[t_] := -1.4*10^3*B15^-2*v3[t]^-2(*(10)式，单位s*)
tacc[t_] := II[t]*\[CapitalOmega][t]/Nacc[t](*(21)式，单位s*)
tbv[t_] := 1/(1.4*10^-10*R115^5*M14^-1*T9^6*
v3[t]^2*(1 + 284.5*R115^4*v3[t]^4*T9^-2*\[Alpha][t]^2 +
3.16*10^4*R115^8*v3[t]^8*T9^-4*\[Alpha][t]^4 +
1.08*10^6*R115^12*v3[t]^12*T9^-6*\[Alpha][t]^6))(*(4)式*)
Aplus[t_] := 1 + (3*\[Alpha][t]^2*Jtilder)/(2*Itilder)(*(28)式下面*)
Aminus[t_] := 1 - (3*\[Alpha][t]^2*Jtilder)/(2*Itilder)(*(28)式下面*)
d\[Alpha][
t_] := \[Alpha][
t]*((Mbdot[t]*Aminus[t])/(2*M[t]*Aplus[t]) -
Aminus[t]/Aplus[t]*(1/tsv + 1/tbv[t] + 1/tBt[t]) -
1/Aplus[t]*(1/tdip[t] + 1/tacc[t] + 1/tBgw[t]) - 1/
tgw[t])(*\[Alpha]'，(28)式*)
d\[CapitalOmega][
t_] := \[CapitalOmega][
t]*(-(Mbdot[t]/(Aplus[t]*M[t])) +
1/Aplus[t]*(1/tdip[t] + 1/tacc[t] + 1/tBgw[t]) - (
3*\[Alpha][t]^2*Jtilder)/(
Itilder*Aplus[t])*(1/tsv + 1/tbv[t] + 1/
tBt[t]))(*\[CapitalOmega]，(27)式*)
dBt[t_] := (4/(3*Pi))^(1/2)*
B*\[Alpha][t]^2*\[CapitalOmega][t](*Subscript[B, t]，(6)式*)
NDSolve[{Mb'[t] == Mbdot[t], \[Alpha]'[t] ==
d\[Alpha][t], \[CapitalOmega]'[t] == d\[CapitalOmega][t],
Bt'[t] == dBt[t],
Mb[0] == Mbi, \[Alpha][0] == 10^-8, \[CapitalOmega][0] == (2*Pi)/(
3*10^3), Bt[0] == 100}, {Mb[t], \[Alpha][t], \[CapitalOmega][t],
Bt[t]}, {t, 0, 10000}]
还是不行啊，前面的问题是遇到无穷大，忽略了也能得到图形，但是加上后面的方程组成微分方程组就完全不行了
