search

SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
556 stars 210 forks source link

IMEX instability with SplitODEProblem #930

Open ChrisRackauckas opened 5 years ago

ChrisRackauckas commented 5 years ago 
using ApproxFun, Sundials, Plots; gr()

S = Fourier()
n = 100
x = points(S, n)
T = ApproxFun.plan_transform(S, n)
Ti = ApproxFun.plan_itransform(S, n)

# Convert the initial condition to Fourier space
u₀ = T*cos.(cos.(x.-0.1))
D2 = Derivative(S,2)
L = D2[1:n,1:n]

# By making it Diagonal, DiffEq internally specializes the linear solver
using DiffEqOperators, LinearAlgebra, DifferentialEquations

A = DiffEqArrayOperator(Diagonal(L))

function f(dû,û,tmp,t)
  # Transform u back to point-space
  mul!(tmp,Ti,û)
  # apply nonlinear function 0.75sqrt(u)-u in point-space
  @. tmp = 0.75sqrt(tmp) - tmp
  mul!(dû,T,tmp) # Transform back to Fourier space
end

# Define u' = Au + f
prob = SplitODEProblem(A, f, u₀, (0.0,10.0), similar(u₀));

sol = solve(prob, KenCarp4())
plot(x,Ti*sol(0.0))
plot!(x,Ti*sol(0.5))
plot!(x,Ti*sol(1.0))
plot!(x,Ti*sol(2.0))
ChrisRackauckas commented 5 years ago

I think this is another instance of the linearity specializations not working out.

ChrisRackauckas commented 5 years ago

@YingboMa you might want to take a look at this one. It is likely related to https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/839 . Somewhere when the Jacobian re-use logic got reduce I think a bug was introduced. These all need to become regression tests.

ChrisRackauckas commented 5 years ago

Capture Looks like the sign is flipped. Maybe the linear solver handling didn't flip the sign in the linear case when it was changed to W_transform style?