Question: DifferentialEquations.jl and Runge-Kutte RK4 error

General Usage
,
1

I’m trying to program a simple matrix-Ricatti ODE solver via the DifferentialEquations.jl package and the RK4 algorithm (I tried the default algorithm as well, to no avail). To benchmark my solver, I’m comparing it to the (user-supplied) MATLAB toolbox “CFH Toolbox”, specifically the function cfaffine.m. The benchmark uses a standard explicit Runge-Kutta 4th-order scheme in a straightforward finite-difference implementation (i.e., no other packages are used). Thus, this would seem like an almost apples-to-apples comparison.

However, in my Julia code I start getting errors for initial values with larger imaginary components, something the MATLAB code with the same algorithm handles without issues. Specifically, the errors are “Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.” These occur regardless of imposing a fixed stepsize (however small or large) or not. Further, even for the intial values that compute without error, as soon as the initial values become complex there is a divergence from the MATLAB benchmark.

What am I misunderstanding about the RK4() implementation? Given that the ODE is numerically straightforward, and especially with a forced fixed stepsize, I don’t understand the divergence from the comparable MATLAB solver or it aborting. Here is my Julia code:

# Julia Code
# Ricatti ODE with Runge-Kutta 4th order
# "dt <= dtmin" error occurs with or without fixed stepsize
using DifferentialEquations
using LinearAlgebra

# Original state-space via Duffie-Pan-Singleton 2000 ECMA
K0	= [.02;.01]
K1	= [-1.0 0.0; 0.0 -1.0]
H0	= [0.0 0.0; 0.0 0.0]
H1	= reshape([.16 0.0 0.0 0.0 0.0 0.0 0.0 .09],(2,2,2)) # CUBE
rho0= .005
rho1= [1.0;1.0]

# Quadratic part in Ricatti ODE below
cube_multiply(cube::Array{Float64,3},u) = dropdims(mapslices(i-> u'*i*u, cube; dims=[1,2]);dims=(1,2))

# Ricatti ODE: Stacked system of alpha and beta
function ODE_alphabeta!(du,u,p,t)
	du[1] = rho0 - K0'*u[2:end]-1/2*u[2:end]'*H0*u[2:end]; # alpha (singleton dimension)
	du[2:end] .= rho1 - K1'*u[2:end] - 1/2*cube_multiply(H1,u[2:end]); # beta
	return du
end

function psi(a,x,T)
	tspan = (T,0.0) # solving backwards from T to 0
	uT    = [0;a] # initial values: [alpha(T)=0;beta(T)=a]
	prob  = ODEProblem(ODE_alphabeta!,uT,tspan)
	sol   = solve(prob,RK4();dt=0.005,adaptive = false) # using 4th-order Runge-Kutta with same stepsize as MATLAB implementation; 
       # edited to reflect ChrisRackauckas comment below on adaptive timestep
	alpha = sol.u[end][1]
	beta  = sol.u[end][2:end]
 	return exp(alpha+beta'*x)
end

# Numerical Example with initial conditions
# Eventual goal: integral \int_0^\infty psi(a0+b0*v*im,x0,T)dv so need solvability for large v
a0 = -[1.0,1.0]
b0 = [1.0,1.0]
x0 = [.02,.01]
T0 = 10

psi(a0,x0,T0) # v=0
psi(a0+5.0*b0*im,x0,T0) # v=5
psi(a0+10.0*b0*im,x0,T0) # v=10
psi(a0+50.0*b0*im,x0,T0) # v=50: this throws a Warning: due to "dt<=dtmin. Aborting."

For completeness, below the MATLAB code for replication:

% Benchmark MATLAB code (using package "CFH Toolbox", function "cfaffine.m")
K0 = [.02,.01]
K1 = [-1 0 0 -1]
H0 = [0 0 0 0]
H1 = [.16 0 0 0 0 0 0 .09]
R0 = 0.005
R1 = [1 1]
a0 = -[1;1]
b0 = [1;1]
x0 = [.02 .01]
T0 = 10
% Note that the MATLAB solver is by default in the imaginary space
cfaffine(-a0*i,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1) % same solution as Julia code above
cfaffine(-a0*i+b0*5,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)
cfaffine(-a0*i+b0*10,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)
cfaffine(-a0*i+b0*50,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)
2

What you’re showing is adaptive time stepping. I assume you meant to do: solve(prob,RK4();dt=0.005,adaptive=false)?

1 Like
3

yes, i meant that with adaptive = false. Nevertheless, i get an output of

NaN + NaN*im

when evaluating

psi(a0+50.0*b0*im,x0,T0)

even after adding adaptive = false in the psi function definition, whereas the above-mentioned MATLAB’s package cfaffine.m function

cfaffine(-a0*i+b0*50,x0,T0,K0,K1,H0,H1,R0,R1,[],[],[],1)

happily computes the finite number of

0.3534 + 0.2271i.
4

If you evaluate your ODE RHS on a random vector, do you get the same output in MATLAB and Julia?

5

I found my error – I was unaware that x’ is not the transpose, but the adjoint. i am new to julia, so please forgive me wasting your time with such beginner mistakes.

thus, my solutions coincided for real x, but diverged for complex ones. just for completeness, this code now perfectly replicates the MATLAB benchmark in my initial post (for random vectors and my chosen ones). it also works with dropping the fixed step-size “dt=0.005” and “adaptive=false” settings, so RK4 works as expected.

# minimum viable example: Ricatti ODE with Runge-Kutta 4th order
# Error occurs with or without fixed stepsize
using DifferentialEquations
using LinearAlgebra

# Original state-space via Duffie-Pan-Singleton 2000 ECMA
K0	= [.02;.01]
K1	= [-1.0 0.0; 0.0 -1.0]
H0	= [0.0 0.0; 0.0 0.0]
H1	= reshape([.16 0.0 0.0 0.0 0.0 0.0 0.0 .09],(2,2,2)) # CUBE
rho0= .005
rho1= [1.0;1.0]

# Quadratic part in Ricatti ODE below
cube_multiply(cube::Array{Float64,3},u) = dropdims(mapslices(i-> transpose(u)*i*u, cube; dims=[1,2]);dims=(1,2))

# Ricatti ODE: Stacked system of alpha and beta
function ODE_alphabeta!(du,u,p,t)
	du[1] = rho0 - K0'*u[2:end]-1/2*transpose(u[2:end])*H0*u[2:end]; # alpha (singleton dimension)
	du[2:end] .= rho1 - K1'*u[2:end] - 1/2*cube_multiply(H1,u[2:end]); # beta
	return du
end

function psi(a,x,T)
	tspan = (T,0.0) # solving backwards from T to 0
	uT    = [0;a] # initial values: [alpha(T)=0;beta(T)=a]
	prob  = ODEProblem(ODE_alphabeta!,uT,tspan)
	sol   = solve(prob,RK4();dt=0.005,adaptive = false) # using Runge-Kutta 4th order algorithm with same stepsize as MATLAB implementation
	alpha = sol.u[end][1]
	beta  = sol.u[end][2:end]
 	return exp(alpha+transpose(beta)*x)
end

# Numerical Example with initial conditions
# Eventual goal: integral \int_0^\infty psi(a0+b0*v*im,x0,T)dv so need solvability for large v
a0 = -[1.0,1.0]
b0 = [1.0,1.0]

# a0 = randn(2)
# b0 = randn(2)

x0 = [.02,.01]
T0 = 10.0

psi(a0,x0,T0) # v=0
psi(a0+5.0*b0*im,x0,T0) # v=5
psi(a0+10.0*b0*im,x0,T0) # v=10
psi(a0+50.0*b0*im,x0,T0) # v=50: this throws a Warning: due to "dt<=dtmin. Aborting."

Again, thank you so much and my apologies for such an amateur mistake.

3 Likes