Suggestion for exploiting separable linear-nonlinear ODEs

[LysandrosAn]

In many cases, systems of ordinary differential equations can be separated into a linear and a nonlinear part. This is often encountered in structural dynamics, where an otherwise ideal, elastic structure is interacting with a nonlinear force at specific degrees of freedom. In these cases, the linear section of the ODE alone would create a dense or banded Jacobian matrix (depending on the system convention), while the nonlinear section alone would create a relatively sparse Jacobian structure.

Can we add a further option to the ODEProblem where the user will be able to provide the analytical expression of the linear part (under the assumption that it is known beforehand) together with the expression or just the sparsity of the nonlinear part? If the algorithm possesses this a priori knowledge, then it has to apply finite differences / automatic differentiation solely to the nonlinearities thereby perhaps accelerating the time-integration of the problem. I have tested this in a software for continuation and there is a significant reduction in calculation time, when there are above 30-40 degrees of freedom of which a few are affected by nonlinearities.

Cheers Lysandros Anastasopoulos

[ChrisRackauckas]

This is described https://docs.sciml.ai/latest/types/split_ode_types/ and the solvers specializing on this form are described in https://docs.sciml.ai/latest/solvers/split_ode_solve/#Semilinear-ODE-1 . Let me know if you need anything else.

[ChrisRackauckas]

For reference, he's some benchmarks which make use of it: https://benchmarks.sciml.ai/html/MOLPDE/allen_cahn_fdm_wpd.html

[LysandrosAn]

Τhanks for the references. I am trying to model the Lorenz system by splitting the ODE:

# Load packages and functions
 using LinearAlgebra
 using Plots
 using MAT
 using JLD
 using DifferentialEquations
 using DiffEqOperators

const σ=10.0;
const β=8/3;
const ρ=28.0;

function lin_Term(x,p,t)  # Define linear part
    A=[-σ +σ 0;
       +ρ -1 0;
        0  0 -β]
    DiffEqArrayOperator(A)
end

function nonlin_Term(x,p,t)  # Define non-linear part
    nonlin_Term=[ 0;
                  -x[1]*x[3];
                  +x[1]*x[2]]
end

# Solve differential equation
tspan = (0.0,0.001);
prob = SplitODEProblem(lin_Term, nonlin_Term, [1.0;1.0;0.0], tspan)
sol = solve(prob);

I must be missing something, because it does not work yet.

[ChrisRackauckas]

The first term should not be a function but the operator itself:

 using LinearAlgebra
 using Plots
 using MAT
 using JLD
 using DifferentialEquations
 using DiffEqOperators

const σ=10.0;
const β=8/3;
const ρ=28.0;

A=[-σ +σ 0;
       +ρ -1 0;
        0  0 -β]
DiffEqArrayOperator(A)

function nonlin_Term(x,p,t)  # Define non-linear part
    nonlin_Term=[ 0;
                  -x[1]*x[3];
                  +x[1]*x[2]]
end

# Solve differential equation
tspan = (0.0,0.001);
prob = SplitODEProblem(A, nonlin_Term, [1.0;1.0;0.0], tspan)
sol = solve(prob);