麻烦大家帮我看看,为什么一直运行,就是不出结果。
mathematica吧
全部回复
仅看楼主
level 5
(*Solve the differential equation with the Runge-Kutta 4th order \
method*)
sol = rungekutt4[radialEq2bc, {rh, rout}, {R0, R1}, 0.01];
就是上面这行代码,运行不出来。下面我将完整的代码放上来,请大家帮我看看究竟是怎么回事呢?我简单的说一下,我这个代码的用途:使用四阶龙格库塔法,数值求解散射波函数,点波源的位置在$rs$处并且在施瓦西黑洞左侧,$rh$为黑洞的视界位置,$rout$是空间无穷远。不求代码,只求思路。麻烦大家了。
Quit
ClearAll["Global`*"]
(*Set the values of rin,rout,rs,omega,Veff,a,and l*)
rh = 2 M;
rout = 100 M;
rs = 6 M;
M = 1;
\[Omega] = 2;
l = 2;
a = 1;
(*Define the Runge-Kutta 4th order method*)
rungekutt4[ufunc_, tspan_, y0_, h_] :=
Block[{t0, tn, n, k1, k2, k3, k4, y, t},
If[Dimensions@璐村惂鐢ㄦ埛_00VPR7V馃惥 == {2}, t0 = tspan[[1]]; tn = tspan[[2]],
Return@"error"];
n = Floor[(tn - t0)/h];
t = ConstantArray[0, n + 1];
y = ConstantArray[0, {2, n + 1}];
t[[1]] = t0;
y[[;; , 1]] = y0;
Do[t[[i + 1]] = t[[i]] + h;
k1 = ufunc[t[[i]], y[[;; , i]]];
k2 = ufunc[t[[i]] + h/2, y[[;; , i]] + h*k1/2];
k3 = ufunc[t[[i]] + h/2, y[[;; , i]] + h*k2/2];
k4 = ufunc[t[[i]] + h, y[[;; , i]] + h*k3];
y[[;; , i + 1]] = y[[;; , i]] + h*(k1 + 2*k2 + 2*k3 + k4)/6, {i, 1,
n}];
Return@{t, y}]
(*Define the effective potential Veff*)
Veff[x_, M_, l_] := (1 - (2 M)/x)*((l (l + 1))/x^2 + (2 M)/(x^3))
(*Define the initial conditions*)
R0 = 0;
R1 = 0;
(*Define the second-order radial differential equation*)
radialEq2[x_, y_, M_,
l_, \[Omega]_] := {y[[2]], -((\[Omega]^2 - Veff[x, M, l]) y[[1]])}
(*Define the boundary conditions*)
bcRh = {R[rh] == 0, R'[rh] == 1};
bcRs = {R'[rs] - a*(-1)^l (l + 1/2) R[rs] == 0, R[rs] == rs};
bcRout = {R[rout] == 0, R'[rout] == 1};
(*Define the function for the differential equation and boundary \
conditions*)
radialEq2bc[x_, y_, M_, l_, \[Omega]_] :=
Join[radialEq2[x, y, M, l, \[Omega]], bcRh, bcRs, bcRout]
(*Solve the differential equation with the Runge-Kutta 4th order \
method*)
sol = rungekutt4[radialEq2bc, {rh, rout}, {R0, R1}, 0.01];
Rsol = sol[[2, All, 1]];
RpSol = sol[[2, All, 2]];
(*Extract the solution and plot it*)
Rsol = sol[[2, 1]]
ListLinePlot[Transpose[{sol[[1]], Rsol}]]
2023年04月15日 00点04分 1
吧务
level 15
radialEq2bc的定义明显不符合rungekutt4里ufunc的语法要求。ufunc必须是2参数的。
2023年05月06日 05点05分 2
1