该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

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}]

还是不行啊，前面的问题是遇到无穷大，忽略了也能得到图形，但是加上后面的方程组成微分方程组就完全不行了