update the functions

Author: ChrisRackauckas

benchmarks/NonStiffODE/NLS_wpd.jmd

@@ -33,7 +33,6 @@ The system consists of two coupled (1+1)D nonlinear PDEs:
 
 ```julia
 using OrdinaryDiffEq, DiffEqDevTools, Plots, StaticArrays, FFTW, LinearAlgebra
-using SciMLBenchmarks
 
 # Parameters
 const β₂F = -1.0
@@ -48,6 +47,7 @@ const s4 = β₂F / abs(β₂S)
 const s5 = γc / γs
 const s6 = γf / γs
 const q = 10.0
+const steps=2000.0
 
 # Grid parameters
 const N = 2048
@@ -59,102 +59,97 @@ const T_max = N / 2
 const dT = 1.0
 const T = collect((-N/2):(N/2-1)) * dT
 const ω = 2π * fftfreq(N, 1/dT)
+const α = 0.5*sign(β₂S).*ω.^2.0
 
-# Initial condition functions
-function soliton(q, t; shift=0.0)
-    return sqrt(2*q) * sech.(sqrt(q) * (t .- shift)) * exp.(1im * (t .- shift))
+function f(x,y,z,a)
+    return a.*exp.(-((x.-y).^2.0)./(2*z^2.0))
 end
 
-function f_pulse(t, t0, A, σ²)
-    return A * exp.(-0.5 * (t .- t0).^2 / σ²)
+function soliton(q,t,shift)
+    sqrt(2.0*q)*sech.(sqrt(2.0*q).*(t.-shift))
 end
 
-# Initial conditions
-ψ₀ = soliton(q, T, shift=20.0)
-F₀ = f_pulse(T, -20.0, 10, 1e-18) .* exp.(1.12157im * T)
+t=(-N/2.0:N/2.0-1.0)*dt 
+shift = 20.0
 
-# Combine into state vector: [Re(ψ), Im(ψ), Re(φ), Im(φ)]
-u₀ = vcat(real(ψ₀), imag(ψ₀), real(F₀), imag(F₀))
+psi_0 = soliton(q,t,shift)
+F_0 = f.(t,-20.0,1e-1,10^-2).*cis.(1.12157 * t)
 
-function nls_rhs!(du, u, p, X)
-    # Extract real and imaginary parts
-    ψ_re = @view u[1:N]
-    ψ_im = @view u[N+1:2N]  
-    φ_re = @view u[2N+1:3N]
-    φ_im = @view u[3N+1:4N]
-    
-    # Reconstruct complex arrays
-    ψ = complex.(ψ_re, ψ_im)
-    φ = complex.(φ_re, φ_im)
-    
-    # FFT to frequency domain
-    ψ_hat = fft(ψ)
-    φ_hat = fft(φ)
-    
-    # Apply dispersion terms in frequency domain
-    # First equation: -i∂ψ/∂X = -(sgn(β₂S)/2)(∂²ψ/∂τ²) + nonlinear terms
-    disp_ψ = -1im * (sign(β₂S)/2) * (1im * ω).^2 .* ψ_hat
-    ψ_disp = ifft(disp_ψ)
-    
-    # Second equation dispersion: -(s₄/2)(∂²φ/∂τ²) + is₁(∂φ/∂τ)  
-    disp_φ = -1im * ((s4/2) * (1im * ω).^2 .* φ_hat + s1 * (1im * ω) .* φ_hat)
-    φ_disp = ifft(disp_φ)
-    
-    # Nonlinear terms
-    ψ_mod² = abs2.(ψ)
-    nonlin_ψ = -1im * ψ_mod² .* ψ
-    nonlin_φ = -1im * (s2 * ψ .* conj.(φ) * exp(1im * (s3 + q) * X) + s5 * ψ_mod² .* φ)
+# Split complex initial conditions into real and imaginary parts
+psi_0_fft = fft(psi_0)
+F_0_fft = fft(F_0)
+ini = [real(psi_0_fft); imag(psi_0_fft); real(F_0_fft); imag(F_0_fft)]
+
+## Solution
+
+function ODE!(dX,X,p,t)
+    s2, s3, s5, s6, alpha, Omega, N = p
+    exp_term_psi = cis.(alpha .* t)
+    exp_term_phi = cis.(Omega .* t) 
+
+    # Reconstruct complex arrays from real and imaginary parts
+    A = complex.(X[1:N], X[N+1:2*N]) .* exp_term_psi
+    B = complex.(X[2*N+1:3*N], X[3*N+1:4*N]) .* exp_term_phi
+    psi = ifft(A)
+    phi = ifft(B)
 
-    # Total RHS for each equation
-    dψ_dt = ψ_disp + nonlin_ψ
-    dφ_dt = φ_disp + nonlin_φ
+    rhs1 = 0.5 * cis.(-s3 * t) * s2 .* (phi.^2.0) .+ psi .* abs2.(psi) .+ psi .* s5 .* abs2.(phi)
+    rhs2 = cis.(s3 * t) * s2 .* psi .* conj(phi) .+ phi * s5 .* abs2.(psi) .+ phi .* s6 .* abs2.(phi)
+    rhs1_fft = fft(rhs1)
+    rhs2_fft = fft(rhs2)
 
-    # Store results in du
-    du[1:N] .= real.(dψ_dt)
-    du[N+1:2N] .= imag.(dψ_dt)
-    du[2N+1:3N] .= real.(dφ_dt)
-    du[3N+1:4N] .= imag.(dφ_dt)
+    # Compute complex derivatives and split back into real/imaginary parts
+    dpsi_dt = 1im .* rhs1_fft ./ exp_term_psi
+    dphi_dt = 1im .* rhs2_fft ./ exp_term_phi
 
-    return nothing
+    dX[1:N] = real.(dpsi_dt)         # Real part of psi derivative
+    dX[N+1:2*N] = imag.(dpsi_dt)     # Imaginary part of psi derivative  
+    dX[2*N+1:3*N] = real.(dphi_dt)   # Real part of phi derivative
+    dX[3*N+1:4*N] = imag.(dphi_dt)   # Imaginary part of phi derivative
 end
 
-# Create the ODE problem
-prob = ODEProblem(nls_rhs!, u₀, (0.0, L))
-probstatic = ODEProblem{false}(nls_rhs!, u₀, (0.0, L))
+p = (s2, s3, s5, s6, α, ω, N)
+prob = ODEProblem(ODE!,ini,(0.0,L),p)
 
 # Generate reference solution
-ref_sol = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12, saveat=0.1)
+ref_sol = solve(prob, Vern7(), abstol=1e-12, reltol=1e-12)
 
 # Store solutions for work-precision testing
-probs = [prob, probstatic]
-test_sol = [ref_sol, ref_sol]
+probs = [prob]
+test_sol = [ref_sol]
 
 abstols = 1.0 ./ 10.0 .^ (5:8)
 reltols = 1.0 ./ 10.0 .^ (2:5)
 ```
 
 ```julia
 # Plot initial conditions
-ψ_init = complex.(u₀[1:N], u₀[N+1:2N])
-φ_init = complex.(u₀[2N+1:3N], u₀[3N+1:4N])
+u₀ = prob.u0
+ψ_init = complex.(u₀[1:N], u₀[N+1:2*N])  # Reconstruct from real/imaginary parts
+φ_init = complex.(u₀[2*N+1:3*N], u₀[3*N+1:4*N])
 
-p1 = plot(T, abs.(ψ_init), label="|ψ₀|", title="Initial Conditions", lw=2)
-plot!(p1, T, abs.(φ_init), label="|φ₀|", lw=2)
+p1 = plot(t, abs.(ifft(ψ_init)), label="|ψ₀|", title="Initial Conditions", lw=2)
 xlabel!(p1, "τ")
 ylabel!(p1, "Amplitude")
-xlim!(p1, (-50, 50))
-
-# Plot final solution
-ψ_final = complex.(ref_sol.u[end][1:N], ref_sol.u[end][N+1:2N])  
-φ_final = complex.(ref_sol.u[end][2N+1:3N], ref_sol.u[end][3N+1:4N])
 
-p2 = plot(T, abs.(ψ_final), label="|ψ(L)|", title="Final Solutions", lw=2)
-plot!(p2, T, abs.(φ_final), label="|φ(L)|", lw=2)
+p2 = plot(t, abs.(ifft(φ_init)), label="|φ₀|", lw=2)
 xlabel!(p2, "τ")
 ylabel!(p2, "Amplitude")
-xlim!(p2, (-50, 50))
 
-plot(p1, p2, layout=(1,2), size=(800, 400))
+# Plot final solution
+u_final = ref_sol.u[end]
+ψ_final = complex.(u_final[1:N], u_final[N+1:2*N])
+φ_final = complex.(u_final[2*N+1:3*N], u_final[3*N+1:4*N])
+
+p3 = plot(t, abs.(ifft(ψ_final)), label="|ψ(L)|", title="Final Solutions", lw=2)
+xlabel!(p3, "τ")
+ylabel!(p3, "Amplitude")
+
+p4 = plot(t, abs.(ifft(φ_final)), label="|φ(L)|", lw=2)
+xlabel!(p4, "τ")
+ylabel!(p4, "Amplitude")
+
+plot(p1, p2, p3, p4, layout=(2,2), size=(800, 400))
 ```
 
 ## Low Order Methods
@@ -165,33 +160,18 @@ setups = [
     Dict(:alg=>BS5()),
     Dict(:alg=>Tsit5()),
     Dict(:alg=>Vern6()),
-    Dict(:alg=>Tsit5(), :prob_choice => 2),
-    Dict(:alg=>Vern6(), :prob_choice => 2)
-]
-
-wp = WorkPrecisionSet(probs, abstols, reltols, setups;
-                     save_everystep=false, appxsol=test_sol, maxiters=Int(1e5),
-                     numruns=10)
-plot(wp, title="NLS System - Low Order Methods")
-```
-
-## Higher Order Methods
-
-```julia
-setups = [
     Dict(:alg=>DP8()),
     Dict(:alg=>Vern7()),
     Dict(:alg=>Vern8()),
     Dict(:alg=>Vern9()),
-    Dict(:alg=>Vern7(), :prob_choice => 2),
-    Dict(:alg=>Vern8(), :prob_choice => 2),
-    Dict(:alg=>Vern9(), :prob_choice => 2)
+    Dict(:alg=>ROCK2()),
+    Dict(:alg=>ROCK4()),
 ]
 
 wp = WorkPrecisionSet(probs, abstols, reltols, setups;
                      save_everystep=false, appxsol=test_sol, maxiters=Int(1e5),
                      numruns=10)
-plot(wp, title="NLS System - Higher Order Methods")
+plot(wp)
 ```
 
 ## Adaptive vs Fixed Step Methods
@@ -211,31 +191,11 @@ wp = WorkPrecisionSet(probs, abstols, reltols, setups;
 plot(wp, title="NLS System - Adaptive vs Fixed Step")
 ```
 
-## Special Methods for Oscillatory Problems
-
-```julia
-abstols = 1.0 ./ 10.0 .^ (6:9)
-reltols = 1.0 ./ 10.0 .^ (3:6)
-
-setups = [
-    Dict(:alg=>Vern7()),
-    Dict(:alg=>Vern8()),
-    Dict(:alg=>Vern9()),
-    Dict(:alg=>DP8()),
-    Dict(:alg=>TsitPap8())
-]
-
-wp = WorkPrecisionSet(probs, abstols, reltols, setups;
-                     save_everystep=false, appxsol=test_sol, maxiters=Int(1e5),
-                     numruns=10)
-plot(wp, title="NLS System - Oscillatory Methods")
-```
-
 ### Conclusion
 
 This benchmark demonstrates the performance of various ODE solvers on a large (8192-dimensional) system derived from a coupled nonlinear Schrödinger equation via spectral discretization. The system exhibits:
 
-- High dimensionality (N=2048 complex variables → 4096 real variables)
+- High dimensionality (N=2048 complex variables → 8192 real variables after real/imaginary splitting)
 - Nonlinear dynamics with quadratic and cubic nonlinearities
 - Oscillatory behavior due to dispersion and phase terms
 - Conservation properties (energy, momentum)

benchmarks/NonStiffODE/Project.toml

@@ -1,8 +1,11 @@
 [deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LSODA = "7f56f5a3-f504-529b-bc02-0b1fe5e64312"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ODEInterface = "54ca160b-1b9f-5127-a996-1867f4bc2a2c"
 ODEInterfaceDiffEq = "09606e27-ecf5-54fc-bb29-004bd9f985bf"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
@@ -12,9 +15,12 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SciMLBenchmarks = "31c91b34-3c75-11e9-0341-95557aab0344"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
+TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
 
 [compat]
+BenchmarkTools = "1"
 DiffEqDevTools = "2.30"
+DifferentialEquations = "7"
 FFTW = "1"
 LSODA = "0.7"
 ODEInterface = "0.5"
@@ -25,3 +31,4 @@ Plots = "1.4"
 SciMLBenchmarks = "0.1"
 StaticArrays = "1"
 Sundials = "4.2"
+TerminalLoggers = "0.1"