ModelingToolkit.jl icon indicating copy to clipboard operation
ModelingToolkit.jl copied to clipboard

Initialization systems don't support parametric guesses

Open MarcBerliner opened this issue 1 year ago • 0 comments

For example,

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t
D = Differential(t)

@parameters a = -1.0
@variables x(t) [guess = a]

eqs = [D(x) ~ a]

initialization_eqs = [1 ~ exp(1 + x)]

@named sys = ODESystem(eqs, t; initialization_eqs)
sys = complete(structural_simplify(sys))

tspan = (0.0, 0.2)
prob = ODEProblem(sys, [], tspan, [])

@show prob.f.initializeprob.u0

sol = solve(prob.f.initializeprob; show_trace=Val(true))

the exact solution to the initialization_eqs is x(0) = -1. If we show the trace, you can see that it begins with an initial value of x(0) = 0

prob.f.initializeprob.u0 = [0.0]

Algorithm: NewtonRaphson(
   descent = NewtonDescent()
)

----     -------------        -----------         
Iter     f(u) inf-norm        Step 2-norm         
----     -------------        -----------         
0        1.71828183e+00       6.46802053e-314     
...
6        0.00000000e+00       1.22346577e-12      
Final    0.00000000e+00      
----------------------      
retcode: Success
u: 1-element Vector{Float64}:
 -1.0

When we change the guess for x(t) to a Float64, i.e., @variables x(t) [guess = -1.0], it respects the guess:

prob.f.initializeprob.u0 = [-1.0]

Algorithm: NewtonRaphson(
   descent = NewtonDescent()
)

----     -------------        -----------         
Iter     f(u) inf-norm        Step 2-norm         
----     -------------        -----------         
0        0.00000000e+00       4.94065646e-324     
1        0.00000000e+00       0.00000000e+00      
Final    0.00000000e+00      
----------------------      
retcode: Success
u: 1-element Vector{Float64}:
 -1.0

cc: @AayushSabharwal

MarcBerliner avatar May 15 '24 15:05 MarcBerliner