It should be something simple, but I really cannot spot the error.
Lets suppose I want to solve a simple mass-spring-damper problem
m D2(x) + c D(x) + k x = f(t)
Following MTK documentation I can use
@parameters t k m ζ
@variables x(t) F(t)
D = Differential(t)
D2 = D*D
f(t) = 100.0/m;
eq = [ F ~ f(t),
D2(x) ~ F - (k/m)*x - ( ζ*2*sqrt(k*m)*D(x) )/m ]
sys = ODESystem(eq,t)
sys_simplified = structural_simplify(sys)
sys_lower = ode_order_lowering(sys_simplified);
so I can keep track of x, F and D(x), as expected. It is working.
Now suppose I want to keep track of the damping force (actually, I want to use it in a more complex operation to define the RHS of D2(x)). I can try
@parameters t k m ζ
@variables x(t) F(t) FD(t)
D = Differential(t)
D2 = D*D
f(t) = 100.0/m
eq = [ F ~ f(t),
FD ~ ( ζ*2*sqrt(k*m)*D(x) )/m
D2(x) ~ F - (k/m)*x - FD ]
sys = ODESystem(eq,t)
sys_simplified = structural_simplify(sys)
sys_lower = ode_order_lowering(sys_simplified);
but when I try to use this second approach I get
prob = ODEProblem(sys_lower,u0,tspan,p,jac=true);
ArgumentError: Any[FD(t)] are missing from the variable map.
I read the documentation and it says that MTK can miss some variables in the mapping, and the suggestion is to use
ODESystem(eq,t).
It is not working in this case, so I suspect it is something related to the presence of D(x) in FD. But at the same time, the first version works and also has D(x) in the RHS.
The complete code to obtain the error is
using LinearAlgebra, Symbolics, ModelingToolkit,
using Plots, DifferentialEquations
@parameters t k m ζ
@variables x(t) F(t) FS(t) FD(t)
D = Differential(t)
D2 = D*D
f(t) = 100.0/m
fd(t) = ( ζ*2*sqrt(k*m)*D(x) )/m
eq = [ F ~ f(t),
FD ~ ( ζ*2*sqrt(k*m)*D(x) )/m,
D2(x) ~ F -(k/m)*x - FD ]
sys = ODESystem(eq,t)
sys_simplified = structural_simplify(sys)
sys_lower = ode_order_lowering(sys_simplified)
u0 = [D(x) => 0.0 ,
x => 1.0 ,
]
p = [ k => 10.0,
m => 1.0,
ζ => 0.1
]
tspan = (0.0, 10.0)
prob = ODEProblem(sys_lower,u0,tspan,p)
sol = solve(prob)
Thank your for your assistance.