Solving system of spatial Ordinary Differential Equations

Specific Domains Modelling & Simulations
, ,
1

Hello all,

I am Federico, and an undergraduate student approaching Julia for the first time as a part of my studies in computational science. I am currently struggling with solving a system of coupled Ordinary Differential Equations dependent on x, y. I have a system of non-linear PDEs, and using MethodOfLines.jl and NonLinearSolve.jl yields a vector of NaNs, with a return code Unstable, as shown below:

retcode: Unstable
Interpolation: Dict{Num, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}}
ivs: 1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 xdomain:(0.0:0.04:1.0,)
u: Dict{Num, Vector{Float64}} with 4 entries:
  B(x) => [0.0, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  NaN, -Inf, -2.2…
  C(x) => [0.0, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  NaN, Inf, -1.23…
  A(x) => [0.0, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  NaN, Inf, -5.53…
  D(x) => [0.0, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  NaN, -Inf, -6.5

However, when solving this system manually, manufacturing a system of discrete equations, with an approximation order of 2 the NonLinearSolve.jl using NewtonRaphson() is solving perfectly, with a retcode Success. Moreover, in the output above, I can read Interpolations.GriddedInterpolation which is making me think the system is getting simplified and thus returning an unstable solution. This is a simplified example showing this problem using MethodOfLines.jl.

using ModelingToolkit
using NonlinearSolve
using MethodOfLines
using Symbolics
using DomainSets

@parameters x, t;

@variables A(..), B(..), C(..), D(..);

# 4 harmonics nonlinear
eqs = Equation[-100A(x) - 50exp(-40(x^2)) - 9Differential(x)(Differential(x)(A(x))) - (12000//1)*(D(x)^2)*B(x) - (1500//1)*(A(x)^2)*B(x) - (1500//1)*(B(x)^3) - (12000//1)*B(x)*(C(x)^2) ~ 0, -100B(x) - 9Differential(x)(Differential(x)(B(x))) + (12000//1)*(D(x)^2)*A(x) + (1500//1)*(A(x)^3) + (1500//1)*A(x)*(B(x)^2) + (12000//1)*A(x)*(C(x)^2) ~ 0, -9Differential(x)(Differential(x)(C(x))) - 400C(x) - (12000//1)*(D(x)^3) - (6000//1)*D(x)*(A(x)^2) - (6000//1)*D(x)*(B(x)^2) - (12000//1)*D(x)*(C(x)^2) ~ 0, -400D(x) - 9Differential(x)(Differential(x)(D(x))) + (12000//1)*(D(x)^2)*C(x) + (6000//1)*(A(x)^2)*C(x) + (6000//1)*(B(x)^2)*C(x) + (12000//1)*(C(x)^3) ~ 0]

bcs = Equation[A(0.0) ~ 0.0, A(1.0) ~ 0.0, B(0.0) ~ 0.0, B(1.0) ~ 0.0, C(0.0) ~ 0.0, C(1.0) ~ 0.0, D(0.0) ~ 0.0, D(1.0) ~ 0.0]

domain_specs = [x ∈ Interval(0.0, 1.0)]


@named pdesys = PDESystem(eqs, bcs, domain_specs, collect([:A, :B, :C, :D]), (A(x), B(x), C(x), D(x)))

discretization = MOLFiniteDifference(x => 0.04, nothing; approx_order=2)

prob = discretize(pdesys, discretization)

sol = solve(prob, NewtonRaphson(); reltol=1e-5, abstol=1e-5)

Is there any way to have MethodOfLines.jl to solve the entire problem, without this simplification?

Thank you for your time and help.