diff --git a/src/perform_step/symplectic_perform_step.jl b/src/perform_step/symplectic_perform_step.jl index 44a16bb110..3e4f6f849f 100644 --- a/src/perform_step/symplectic_perform_step.jl +++ b/src/perform_step/symplectic_perform_step.jl @@ -568,7 +568,7 @@ end kdu = f.f1(du, u, p, tnew) du = du + dt * a5 * kdu - tnew = tnew + t + a5 * dt + tnew = tnew + a5 * dt ku = f.f2(du, u, p, tnew) u = u + dt * b6 * ku @@ -625,7 +625,7 @@ end f.f1(kdu, du, u, p, tnew) @.. broadcast=false du=du + dt * a5 * kdu - tnew = tnew + t + a5 * dt + tnew = tnew + a5 * dt f.f2(ku, du, u, p, tnew) @.. broadcast=false u=u + dt * b6 * ku diff --git a/test/algconvergence/symplectic_tests.jl b/test/algconvergence/symplectic_tests.jl index 3d496ee4c8..c3cf9145b9 100644 --- a/test/algconvergence/symplectic_tests.jl +++ b/test/algconvergence/symplectic_tests.jl @@ -67,3 +67,29 @@ end @test_throws ArgumentError solve(prob, alg(); dt = dt) end end + +function motionfuncDirect1(dv,v,u,p,t) + # 1:Electron, 2: Be + ω_1,ω_2,γ,m_1,m_2,η,ω_d=p + dv[1]=-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1 + dv[2]=-ω_2^2*u[2]-γ*u[1]/m_2 +end + +function motionfuncDirect1(v,u,p,t) + # 1:Electron, 2: Be + ω_1,ω_2,γ,m_1,m_2,η,ω_d=p + [-ω_1^2*u[1]*(1+η*cos(ω_d*t))-γ*u[2]/m_1,-ω_2^2*u[2]-γ*u[1]/m_2] +end + +param=[90386.15717208837, 3938.9288690708827, 8560.718748264337, 0.000544617021484666, 8.947079933513658, 0.7596480420227258, 78778.57738141765] +u0_direct=zeros(2) # mm, mm +v0_direct = [0.0, 135.83668926684385] +tspan=(0.0, 1.321179076090661) +prob_direct=SecondOrderODEProblem(motionfuncDirect1,v0_direct,u0_direct,tspan,param) +dt=2e-8 +ref=solve(prob_direct,DPRKN12(),abstol=1e-12,reltol=1e-12,maxiters=1e7,saveat=0.01) + +@testset "symplectic time-dependent $alg-$iip-$pa" for (alg, x, d) in ALGOS + sol=solve(prob_direct,alg,dt=dt,saveat=0.01) + @test maximum(ref-sol) < 1e-3 +end