RK4 -> Tsit5
This commit is contained in:
parent
95c1083693
commit
63feefbd1b
|
|
@ -27,19 +27,19 @@ end
|
||||||
divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
|
divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
|
||||||
shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
|
shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
|
||||||
the other parameters are the same from dirac!(...)."
|
the other parameters are the same from dirac!(...)."
|
||||||
function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true)
|
function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, r_max, divs; shooting=true, normalize=true, algo=Tsit5())
|
||||||
Δr = r_max / divs
|
Δr = r_max / divs
|
||||||
|
|
||||||
if shooting
|
if shooting
|
||||||
@assert divs % 2 == 0 "divs must be an even number when shooting=true"
|
@assert divs % 2 == 0 "divs must be an even number when shooting=true"
|
||||||
prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2))
|
prob = ODEProblem(dirac!, [0, 1], (r_max, r_max / 2))
|
||||||
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
|
sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
|
||||||
wf_right = reverse(hcat(sol.u...); dims=2)
|
wf_right = reverse(hcat(sol.u...); dims=2)
|
||||||
r_max = r_max / 2 # for the next step
|
r_max = r_max / 2 # for the next step
|
||||||
end
|
end
|
||||||
|
|
||||||
prob = ODEProblem(dirac!, [0, 1], (0, r_max))
|
prob = ODEProblem(dirac!, [0, 1], (0, r_max))
|
||||||
sol = solve(prob, RK4(), p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
|
sol = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=Δr)
|
||||||
wf = hcat(sol.u...)
|
wf = hcat(sol.u...)
|
||||||
|
|
||||||
if shooting # join two segments
|
if shooting # join two segments
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue