Linear, non-adaptive time-dependent PDE discretization runs into Newton iteration

[dkarrasch]

As already discussed a while ago on Gitter, we encountered some problems solving linear time-dependent ODEs. MWE:

using OrdinaryDiffEq, DiffEqOperators, LinearAlgebra, SparseArrays
u0 = sin.(range(0, 2pi, length = 10))
tspan = (0.0, 1.0)
f = t-> fill(-100t , 10)
update_coeffs! = (A,u,p,t) ->  (@show t; vals = A.nzval; vals .= f(t); A)
A = DiffEqArrayOperator(spdiagm(0 => fill(0.0, 10)); update_func = update_coeffs!)
prob = ODEProblem(A, u0, tspan)
solve(prob, ImplicitEuler(), adaptive = false, dt = 0.1)

yields

t = 0.0
t = 0.0
t = 0.1
t = 0.1
┌ Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt.
└ @ DiffEqBase ~/.julia/packages/DiffEqBase/6ewdP/src/integrator_interface.jl:168
retcode: ConvergenceFailure
Interpolation: 3rd order Hermite
t: 1-element Array{Float64,1}:
 0.0
u: 1-element Array{Array{Float64,1},1}:
 [0.0, 0.642788, 0.984808, 0.866025, 0.34202, -0.34202, -0.866025, -0.984808, -0.642788, -2.44929e-16]

For non-adaptive ImplicitEuler with a constant stepsize of dt = 0.1, one would expect that nlsolve just solves a linear system of the form (I-dt*A(t))*u = b? In more complicated cases, we saw like more than 50 updates for the same t, probably part of the linearization code, thinking that the DiffEqArrayOperator encodes a quasilinear PDE operator? But we are dealing with truly linear, yet time-dependent operators. Also, the error message disappears if one reduces the -100 to -10 in the definition of f, or reduces the time step to 0.01.

cc: @Phiro333

[ChrisRackauckas]

@YingboMa can you double check the update coefficients logic?

[dkarrasch]

Is there any update on this issue?

[ChrisRackauckas]

There have been a number of changes to the nonlinear solver to prepare for this, but the last step hasn't been done yet. It shouldn't be too hard at this point though. It's whenever I, @devmotion, or @YingboMa have time really. It's the same issue as #914

[dkarrasch]

Cool, great to know there's been progress!

[hersle]

What is the state of exploiting ODE linearity in stiff solvers now in 2025? I am very interested in getting this to work.

[ChrisRackauckas]

First check, is there anything to be done? Newton's method actually has a guarantee that if f is linear then it will step to the root in 1 step regardless of the initial guess. So I would first check whether this property is already had. I believe the way we have the adaptivity algorithm setup, if you converge in one step then it shouldn't ever compute a new jacobian or factorization, so this would be automatically handled since like 2021.

[ChrisRackauckas]

We should make a test using sol.stats and check the number of factorizations and new Jacobians.

[hersle]

Newton's method actually has a guarantee that if f is linear then it will step to the root in 1 step regardless of the initial guess.

Ok, that's great, then it should not be necessary to special-case linear ODEs.

Here is an updated and simplified version of the original example that runs today:

using OrdinaryDiffEq
f(du, u, p, t) = (du .= -100t .* u)
u0 = sin.(range(0, 2π, length = 10))
tspan = (0.0, 1.0)
prob = ODEProblem(f, u0, tspan)

Non-adaptive still fails and it looks like Newton's method does not converge in 1 iteration

julia> sol = solve(prob, ImplicitEuler(nlsolve = NLNewton()); adaptive = false, dt = 0.1)
┌ Warning: Newton steps could not converge and algorithm is not adaptive. Use a lower dt.
└ @ SciMLBase ~/.julia/packages/SciMLBase/YE7xF/src/integrator_interface.jl:671
retcode: ConvergenceFailure
Interpolation: 3rd order Hermite
t: 1-element Vector{Float64}:
 0.0
u: 1-element Vector{Vector{Float64}}:
 [0.0, 0.6427876096865393, 0.984807753012208, 0.8660254037844387, 0.3420201433256689, -0.34202014332566866, -0.8660254037844385, -0.9848077530122081, -0.6427876096865396, -2.4492935982947064e-16]

Adaptive works, but how are the stats here @ChrisRackauckas? There are many nonlinear solver iterations.

julia> sol = solve(prob, ImplicitEuler(nlsolve = NLNewton()); adaptive = true); sol.stats
SciMLBase.DEStats
Number of function 1 evaluations:                  1664
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    368
Number of linear solves:                           1274
Number of Jacobians created:                       2
Number of nonlinear solver iterations:             1264
Number of nonlinear solver convergence failures:   0
Number of fixed-point solver iterations:           0
Number of fixed-point solver convergence failures: 0
Number of rootfind condition calls:                0
Number of accepted steps:                          250
Number of rejected steps:                          117

Restricting Newton to 1 iteration is very inefficient:

julia> sol = solve(prob, ImplicitEuler(nlsolve = NLNewton(max_iter = 1)); adaptive = true); sol.stats
SciMLBase.DEStats
Number of function 1 evaluations:                  710401
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    227519
Number of linear solves:                           227549
Number of Jacobians created:                       43895
Number of nonlinear solver iterations:             183655
Number of nonlinear solver convergence failures:   139756
Number of fixed-point solver iterations:           0
Number of fixed-point solver convergence failures: 0
Number of rootfind condition calls:                0
Number of accepted steps:                          43899
Number of rejected steps:                          0

[ChrisRackauckas]

Number of Jacobians created:                       2

That's 1 more than I'd expect, but close. I'd want to dig into why.

Note @oscardssmith is finishing it up so that NonlinearSolveAlg is used as the default nonlinear solver in OrdinaryDiffEq.jl (#2727). When that lands, I'd hope it just matches performance of the built-in one, and then in v7 we just remove it. So that means... this one is about to be removed, so I'd look to do this optimization on there. The reuse logic is reused though, but it may be easier to do it direct than to rebase it. But it would be good to figure out why more linear solvers.

Note that W = I - gammaJ, where gamma = constdt. So as dt changes you do need to refactorize if adaptive. That's fundamental. So that mostly looks right, other than the fact that it does 2 instead of 1 of J which is ??? maybe it's counting the zero initialization or something.

So that's adaptive. For non-adaptive, not converging is expected on any sufficiently hard problem. The reason is that if you Newton guess diverges, the general strategy is to decrease dt and try again, since then the guess of u_{n+1} = u_n converges to be true, and thus you're guaranteed to get a convergent newton step at some point. However, a bad or unstable newton guess can cause even a nonlinear case to not converge well. That might just be because of the specifics of the ODE-specific Newton variation. It's likely a lot more robust with NonlinearSolveAlg, where certain details of ensuring it's factorizable around a singularity etc. are just handled by NonlinearSolve.jl directly. So, I'd test that again with NonlinearSolveAlg. We should make sure non-adaptive with nonlinearsolvealg, particularly with globalization strategies like trust region, are well-behaved as that will soon be the right way to non-adaptive.

[ChrisRackauckas]

And for the record, finishing this change to NonlinearSolveAlg parity, i.e. #2727 and other steps, is something that really needs some help, so if you're looking to pitch in with nonlinear solver work then that's the right place to be. We really just need to profile it against the built-in one and figure out all of the remaining overheads for why NonlinearSolve.jl is slower. Once it's not slower, then we make it the default. We're close.

[hersle]

Thank you, I will look out for #2727.

For adaptive: the example above has 250+117 (attempted) steps, but 1264 nonlinear solver iterations. Should not 1 attempted step require 1 nonlinear solve which solves in exactly 1 linear solve if the ODE is linear? Or does this not hold because the Jacobian is reused over time steps?

For non-adaptive, not converging is expected on any sufficiently hard problem.

Does not this contradict the property that Newton should converge in 1 iteration on a linear ODE regardless of the initial guess? And if that is not the case, why would one not simply do one direct linear solve in this case to get rid of the problem?