fortran中的射击方法（中子星振荡）

本文介绍了fortran中的射击方法（中子星振荡）的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述 我一直在用fortran 90写剧本，用射击方法解决中子星的径向振荡问题。但由于不明原因，我的程序从未解决。如果没有拍摄方法组件，程序会顺利运行，因为它成功构建了恒星。但是一旦射门进入，一切都会消失。 PROGRAM ROSCILLATION2 USE eos_parameters IMPLICIT NONE INTEGER :: i，j，k ，l INTEGER，PARAMETER :: N_ode = 5 REAL，DIMENSION（N_OD）:: y REAL（8）:: rho0_cgs，rho0，P0，r0，phi0，pi REAL（8）:: r，rend，mass，P，phi，delta，xi，eta REAL（8）:: step，omega，omegastep，tiny，rho_print，Radius，B，a2，s0 ，lamda，E0，E EXTERNAL :: fcn !!!!用户输入 rho0_cgs = 2.D + 15！以cgs为单位的中央密度 step = 1.D-4！步长dr omegastep = 1.D-2！步长d（Ω） tiny = 1.D-8！小数字P（R）/ P（0）来定义星表面 !!!!!!!!! open（unit = 15，file =data.dat，status =new） pi = ACOS（-1.D0） a2 =（（（ （1.6022D-13）** 4）*（6.674D-11）*（（2.997D8）** - 7）*（（1.0546D-34）** - 3）*（1.D6））** （0.5D0））* a2_MeV！转换为编码单位（km ^ -1） B =（（1.6022D-13）** 4）*（6.674D-11）*（（2.997 D8）** - 7）*（（1.0546D-34）** - 3）*（1.D6）* B_MeV！转换为代码单位（km ^ -2） s0 =（ 1.D0 / 3.D0） - （1 /（6 * pi ** 2））* a2 *（（1 /（16 * pi ** 2）* a2 ** 2 +（pi ** - 2）* a1 *（rho0-B））** - 0.5）！在r = 0时声速的平方 lamda = -0.5D0 * log（1-2 * y（1）/ r ） E0 =（r0 ** - 2）* s0 * exp（lamda + 3 * phi0） rho0 = rho0_cgs * 6.67 D-18 / 9.D0！将rho0转换为代码单位（km ^ -2） ！计算中心压力P0 P0 =（1.D0 / 3.D0）* rho0 - （4.D0 / 3.D0）* B - （1.D0 /（a4 *（12.D0 ）*（pi ** 2）））* a2 ** 2 - & &（a2 /（（3.D0）* a4））*（（（1.D0 /（16.D0 * pi ** 4））* a2 ** 2 +（1.D0 /（ pi ** 2））* a4 *（rho0-B））** 0.5D0） !!度量函数phi phi0 = 0.1D0的初始值！任意（需要稍后调整） r0 = 1.D-30！整合起点 ！设置初始条件 !!!!!!!!!!!!!!!!! !!开始整合循环 !!!!!!!!!!!!!!!!! （3）= 1 /（3 * E0）（2）= P0 $ b $由（3）= phi0 $ b $（1）= 0.D0 $ b $ y（5）= 1 omega = 2 * pi * 1000 /（2.997D5）！1kHz以代码单位 DO l = 1,100 omega = omega + omegastep！拍摄方法部分 DO i = 1，1000000000 rend = r0 + REAL（i）* step call oderk（r，rend，y，N_ode，fcn） （4） $ br = rend mass = y（1） P = y（2） phi = y（3） xi = y eta = y（5） IF（P WRITE（*，*）Central density（10 ^ 14 cgs）= rho0_cgs / 1.D14 WRITE（*，*）质量（太阳质量）=，质量/1.477D0 WRITE（*，*）半径（km）=，r WRITE（*，*）紧凑性M / R，质量/ r WRITE（15，*）（Ω* 2.997D5 /（2 * pi）），y（5） GOTO 21 ENDIF ENDDO ENDDO 21继续 结束程序roscillation2 ！ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! $ b $ SUBROUTINE fcn（r，y，yprime） USE eos_parameters IMPLICIT NONE REAL（8），DIMENSION（5）:: y，yprime REAL（8）:: r，m，P，phi，rho，pi，B，a2，xi，eta，W，Q，E，s，lamda，omega INTEGER :: j $（b）= ACOS（-1.D0） a2 =（（（（1.6022D-13）** 4）*（6.674D-11）*（（2.997D8）** - 7） *（（1.0546D-34）** - 3）*（1.D6））**（0.5D0））* a2_MeV！转换为代码单位（km ^ -1） B =（（1.6022D -13）** 4）*（6.674D-11）*（（2.997D8）** -7）*（（1.0546D-34）** -3）*（1.D6）* B_MeV！转换为编码单位（km ^ -2）m = y（1） P = y（2） phi = y（3） xi = y（4） eta = y（5） rho = 3.D0 * P + 4.D0 * B +（（3.D0）/（4.D0 * a4 *（pi ** 2））） * A2 ** 2 +（A2 / A4）*&安培; （（3.D0 /（pi ** 2））* a4 * b2 +（（（（D.D0）*（pi ** 4）））* a2 ** 2 + （P + B）））** 0.5D0） s =（1.D0 / 3.D0） - （1 /（6 * pi ** 2））* a2 *（（1 / （16 * pi ** 2）* a2 ** 2 +（pi ** - 2）* a4 *（rho-B））** - 0.5）！声速的平方 W =（r ** - 2）*（rho + P）* exp（3 * lamda + phi） E =（r ** - 2）* s * exp（lamda + 3 * phi ）$（$） Q =（r ** - 2）* exp（lamda + 3 * phi）*（rho + P）*（（yprime（3）** 2）+ 4 *（r * * -1）* yprime（3）-8 * pi * P * exp（2 * lamda）） yprime（1）= 4.D0 * pi * rho * r ** 2 $ b $ yprime（2）= - （rho + P）*（m + 4.D0 * pi * P * r ** 3）/（r *（r-2.D0 * m）） yprime（3）=（m + 4.D0 * pi * P * r ** 3）/（r *（r-2.D0 * m）） yprime （4）= y（5）/（3 * E） yprime（5）= - （W * omega ** 2 + Q）* y（4） END SUBROUTINE fcn !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ！ ！ Runge-Kutta方法（来自Numerical Recipes）！ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 子程序oderk（ri，re，y，n，derivs） INTEGER，PARAMETER :: NMAX = 16 REAL（8）:: ri，re，step REAL（8），DIMENSION（NMAX）:: y，dydx，yout EXTERNAL :: derivs，rk4 调用derivs（ri，y，dydx） step = re-ri CALL rk4（y，dydx，n，ri，step，yout，derivs） do i = 1，n y（i）= yout（i） enddo 返回结束子程序oderk SUBROUTINE RK4（Y，DYDX，N，X，H，YOUT，DERIVS） INTEGER，PARAMETER :: NMAX = 16 REAL（8）:: H，HH，XH，X，H6 REAL（8），DIMENSION（N）:: Y，DYDX，YOUT REAL（8），DIMENSION （NMAX）:: YT，DYT，DYM EXTERNAL :: derivs HH = H * 0.5D0 H6 = H / 6D0 XH = X + HH DO I = 1，N YT（I）= Y（I）+ HH * DYDX（I） ENDDO CALL DERIVS（XH，YT，DYT） DO I = 1，N YT（I）= Y（I）+ HH * DYT（I） ENDDO CALL DERIVS（XH，YT，DYM） DO I = 1，N YT（I）= Y（I）+ H * DYM（I） DYM（I）= DYT（I）+ DYM（I） ENDDO CALL DERIVS（X + H，YT，DYT） 请问我= 1，N YOUT（I）= Y（I）+ H6 *（DYDX（I）+ DYT（I）+ 2 * DYM（I）） ENDDO END SUBROUTINE RK4 任何答复都会很棒我真的很郁闷调试。解决方案 你的程序因为这一行而炸毁： yprime（5）= - （W * omega ** 2 + Q）* y（4） 子程序 fcn 中的 。在这个子程序中， omega 完全独立于主程序中声明的那个。如果你的编译器足够好（或者被告知）初始化变量，那么这个未初始化并用于表达式，该表达式将包含随机值或零。 如果您希望主程序中的变量 omega 与您在 fcn ，那么你需要以某种方式将该变量传递给 fcn 。由于你设计这个程序的方式，传递它将需要修改你的所有程序来传递 omega ，以便它可以被提供给你所有的 DERIVS （这是您与 fcn 关联的虚拟参数）。 另一种方法是将 omega 放入模块中，使用该模块，您需要访问 omega ，例如在 eos_parameters 中声明它，而不是在 fcn 的作用域单元和你的主程序中声明它。 I have been writing a script in fortran 90 for solving the radial oscillation problem of a neutron star with the use of shooting method. But for unknown reason, my program never works out. Without the shooting method component, the program runs smoothly as it successfully constructed the star. But once the shooting comes in, everything dies. PROGRAM ROSCILLATION2USE eos_parametersIMPLICIT NONEINTEGER ::i, j, k, lINTEGER, PARAMETER :: N_ode = 5REAL, DIMENSION(N_ode) :: yREAL(8) :: rho0_cgs, rho0, P0, r0, phi0, piREAL(8) :: r, rend, mass, P, phi, delta, xi, etaREAL(8) :: step, omega, omegastep, tiny, rho_print, Radius, B, a2, s0, lamda, E0, EEXTERNAL :: fcn!!!! User input rho0_cgs = 2.D+15 !central density in cgs unitstep = 1.D-4 ! step size dromegastep = 1.D-2 ! step size d(omega)tiny = 1.D-8 ! small number P(R)/P(0) to define star surface!!!!!!!!!open(unit=15, file="data.dat", status="new")pi = ACOS(-1.D0)a2 =((((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6))**(0.5D0))*a2_MeV !convert to code unit (km^-1)B = ((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6)*B_MeV !convert to code unit (km^-2)s0 = (1.D0/3.D0) - (1/(6*pi**2))*a2*((1/(16*pi**2)*a2**2 + (pi**-2)*a4*(rho0 - B))**-0.5) !square of the spped of sound at r=0lamda = -0.5D0*log(1-2*y(1)/r)E0 = (r0**-2)*s0*exp(lamda + 3*phi0)rho0 = rho0_cgs*6.67D-18 / 9.D0 !convert rho0 to code unit (km^-2)!! Calculate central pressure P0P0 = (1.D0/3.D0)*rho0 - (4.D0/3.D0)*B - (1.D0/(a4*(12.D0)*(pi**2)))*a2**2 - &&(a2/((3.D0)*a4))*(((1.D0/(16.D0*pi**4))*a2**2+(1.D0/(pi**2))*a4*(rho0-B))**0.5D0)!! initial value for metric function phi phi0 = 0.1D0 ! arbitrary (needed to be adjusted later)r0 = 1.D-30 ! integration starting point!! Set initial conditions!!!!!!!!!!!!!!!!!!!Start integration loop!!!!!!!!!!!!!!!!!r = r0y(1) = 0.D0y(2) = P0y(3) = phi0y(4) = 1/(3*E0)y(5) = 1omega = 2*pi*1000/(2.997D5) !omega of 1kHz in code unit DO l = 1, 1000 omega = omega + omegastep !shooting method partDO i = 1, 1000000000 rend = r0 + REAL(i)*step call oderk(r,rend,y,N_ode,fcn) r = rend mass = y(1) P = y(2) phi = y(3) xi = y(4) eta = y(5) IF (P < tiny*P0) THEN WRITE(*,*) "Central density (10^14 cgs) = ", rho0_cgs/1.D14 WRITE(*,*) " Mass (solar mass) = ", mass/1.477D0 WRITE(*,*) " Radius (km) = ", r WRITE(*,*) " Compactness M/R ", mass/r WRITE(15,*) (omega*2.997D5/(2*pi)), y(5) GOTO 21 ENDIFENDDOENDDO21 CONTINUEEND PROGRAM roscillation2!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!SUBROUTINE fcn(r,y,yprime)USE eos_parametersIMPLICIT NONE REAL(8), DIMENSION(5) :: y, yprimeREAL(8) :: r, m, P, phi, rho, pi, B, a2, xi, eta, W, Q, E, s, lamda, omegaINTEGER :: j pi = ACOS(-1.D0)a2 =((((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6))**(0.5D0))*a2_MeV !convert to code unit (km^-1)B = ((1.6022D-13)**4)*(6.674D-11)*((2.997D8)**-7)*((1.0546D-34)**-3)*(1.D6)*B_MeV !convert to code unit (km^-2)m = y(1)P = y(2)phi = y(3)xi = y(4)eta = y(5)rho = 3.D0*P + 4.D0*B +((3.D0)/(4.D0*a4*(pi**2)))*a2**2+(a2/a4)*&&(((9.D0/((16.D0)*(pi**4)))*a2**2+((3.D0/(pi**2))*a4*(P+B)))**0.5D0)s = (1.D0/3.D0) - (1/(6*pi**2))*a2*((1/(16*pi**2)*a2**2 + (pi**-2)*a4*(rho - B))**-0.5) !square of speed of soundW = (r**-2)*(rho + P)*exp(3*lamda + phi)E = (r**-2)*s*exp(lamda + 3*phi)Q = (r**-2)*exp(lamda + 3*phi)*(rho + P)*((yprime(3)**2) + 4*(r**-1)*yprime(3)- 8*pi*P*exp(2*lamda))yprime(1) = 4.D0*pi*rho*r**2yprime(2) = - (rho + P)*(m + 4.D0*pi*P*r**3)/(r*(r-2.D0*m))yprime(3) = (m + 4.D0*pi*P*r**3)/(r*(r-2.D0*m)) yprime(4) = y(5)/(3*E)yprime(5) = -(W*omega**2 + Q)*y(4)END SUBROUTINE fcn!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Runge-Kutta method (from Numerical Recipes)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!subroutine oderk(ri,re,y,n,derivs) INTEGER, PARAMETER :: NMAX=16REAL(8) :: ri, re, stepREAL(8), DIMENSION(NMAX) :: y, dydx, youtEXTERNAL :: derivs,rk4 call derivs(ri,y,dydx) step=re-ri CALL rk4(y,dydx,n,ri,step,yout,derivs) do i=1,n y(i)=yout(i) enddoreturn end subroutine oderkSUBROUTINE RK4(Y,DYDX,N,X,H,YOUT,DERIVS) INTEGER, PARAMETER :: NMAX=16 REAL(8) :: H,HH,XH,X,H6 REAL(8), DIMENSION(N) :: Y, DYDX, YOUT REAL(8), DIMENSION(NMAX) :: YT, DYT, DYM EXTERNAL :: derivsHH=H*0.5D0 H6=H/6D0 XH=X+HHDO I=1,N YT(I)=Y(I)+HH*DYDX(I) ENDDOCALL DERIVS(XH,YT,DYT) DO I=1,N YT(I)=Y(I)+HH*DYT(I) ENDDOCALL DERIVS(XH,YT,DYM) DO I=1,N YT(I)=Y(I)+H*DYM(I) DYM(I)=DYT(I)+DYM(I) ENDDOCALL DERIVS(X+H,YT,DYT) DO I=1,N YOUT(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2*DYM(I)) ENDDOEND SUBROUTINE RK4Any reply would be great i am just really depressed for the long debugging. 解决方案 Your program is blowing up because of this line:yprime(5) = -(W*omega**2 + Q)*y(4)in subroutine fcn. In this subroutine, omega is completely independent of the one declared in your main program. This one is uninitialized and used in an expression, which will either contain random values or zero, if your compiler is nice enough (or told) to initialize variables. If you want the variable omega from your main program to be the same variable you use in fcn then you need to pass that variable to fcn somehow. Due to the way you've architected this program, passing it would require modifying all of your procedures to pass omega so that it can be provided to all of your calls to DERIVS (which is the dummy argument you are associating with fcn). An alternative would be to put omega into a module and use that module where you need access to omega, e.g. declare it in eos_parameters instead of declaring it in the scoping units of fcn and your main program. 这篇关于fortran中的射击方法（中子星振荡）的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！