Skip to content

Commit ff51e63

Browse files
Merge pull request #2096 from SciML/mcate5
Fix McAte5 time dependence
2 parents c9db428 + f5e0454 commit ff51e63

File tree

2 files changed

+28
-2
lines changed

2 files changed

+28
-2
lines changed

src/perform_step/symplectic_perform_step.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -568,7 +568,7 @@ end
568568
kdu = f.f1(du, u, p, tnew)
569569
du = du + dt * a5 * kdu
570570

571-
tnew = tnew + t + a5 * dt
571+
tnew = tnew + a5 * dt
572572
ku = f.f2(du, u, p, tnew)
573573
u = u + dt * b6 * ku
574574

@@ -625,7 +625,7 @@ end
625625
f.f1(kdu, du, u, p, tnew)
626626
@.. broadcast=false du=du + dt * a5 * kdu
627627

628-
tnew = tnew + t + a5 * dt
628+
tnew = tnew + a5 * dt
629629
f.f2(ku, du, u, p, tnew)
630630
@.. broadcast=false u=u + dt * b6 * ku
631631

test/algconvergence/symplectic_tests.jl

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -67,3 +67,29 @@ end
6767
@test_throws ArgumentError solve(prob, alg(); dt = dt)
6868
end
6969
end
70+
71+
function motionfuncDirect1(dv,v,u,p,t)
72+
# 1:Electron, 2: Be
73+
ω_1,ω_2,γ,m_1,m_2,η,ω_d=p
74+
dv[1]=-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1
75+
dv[2]=-ω_2^2*u[2]-γ*u[1]/m_2
76+
end
77+
78+
function motionfuncDirect1(v,u,p,t)
79+
# 1:Electron, 2: Be
80+
ω_1,ω_2,γ,m_1,m_2,η,ω_d=p
81+
[-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1,-ω_2^2*u[2]-γ*u[1]/m_2]
82+
end
83+
84+
param=[90386.15717208837, 3938.9288690708827, 8560.718748264337, 0.000544617021484666, 8.947079933513658, 0.7596480420227258, 78778.57738141765]
85+
u0_direct=zeros(2) # mm, mm
86+
v0_direct = [0.0, 135.83668926684385]
87+
tspan=(0.0, 1.321179076090661)
88+
prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,param)
89+
dt=2e-8
90+
ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01)
91+
92+
@testset "symplectic time-dependent $alg-$iip-$pa" for (alg, x, d) in ALGOS
93+
sol=solve(prob_direct,alg,dt=dt,saveat=0.01)
94+
@test maximum(ref-sol) < 1e-3
95+
end

0 commit comments

Comments
 (0)