@@ -1552,3 +1552,60 @@ end
1552
1552
expected_tstops = unique! (sort! (vcat (0.0 : 0.075 : 10.0 , 0.1 , 0.2 , 0.65 , 0.35 , 0.45 )))
1553
1553
@test all (x -> any (isapprox (x, atol = 1e-6 ), sol2. t), expected_tstops)
1554
1554
end
1555
+
1556
+ @testset " dae_order_lowering basic test" begin
1557
+ @parameters a
1558
+ @variables x (t) y (t) z (t)
1559
+ @named dae_sys = ODESystem ([
1560
+ D (x) ~ y,
1561
+ 0 ~ x + z,
1562
+ 0 ~ x - y + z
1563
+ ], t, [z, y, x], [])
1564
+
1565
+ lowered_dae_sys = dae_order_lowering (dae_sys)
1566
+ @variables x1 (t) y1 (t) z1 (t)
1567
+ expected_eqs = [
1568
+ 0 ~ x + z,
1569
+ 0 ~ x - y + z,
1570
+ Differential (t)(x) ~ y
1571
+ ]
1572
+ lowered_eqs = equations (lowered_dae_sys)
1573
+ sorted_lowered_eqs = sort (lowered_eqs, by= string)
1574
+ sorted_expected_eqs = sort (expected_eqs, by= string)
1575
+ @test sorted_lowered_eqs == sorted_expected_eqs
1576
+
1577
+ expected_vars = Set ([z, y, x])
1578
+ lowered_vars = Set (unknowns (lowered_dae_sys))
1579
+ @test lowered_vars == expected_vars
1580
+ end
1581
+
1582
+ @testset " dae_order_lowering test with structural_simplify" begin
1583
+ @variables x (t) y (t) z (t)
1584
+ @parameters M b k
1585
+ eqs = [
1586
+ D (D (x)) ~ - b / M * D (x) - k / M * x,
1587
+ 0 ~ y - D (x),
1588
+ 0 ~ z - x
1589
+ ]
1590
+ ps = [M, b, k]
1591
+ default_u0 = [
1592
+ D (x) => 0.0 , x => 10.0 , y => 0.0 , z => 10.0
1593
+ ]
1594
+ default_p = [M => 1.0 , b => 1.0 , k => 1.0 ]
1595
+ @named dae_sys = ODESystem (eqs, t, [x, y, z], ps; defaults = [default_u0; default_p])
1596
+
1597
+ simplified_dae_sys = structural_simplify (dae_sys)
1598
+
1599
+ lowered_dae_sys = dae_order_lowering (simplified_dae_sys)
1600
+ lowered_dae_sys = complete (lowered_dae_sys)
1601
+
1602
+ tspan = (0.0 , 10.0 )
1603
+ prob = ODEProblem (lowered_dae_sys, nothing , tspan)
1604
+ sol = solve (prob, Tsit5 ())
1605
+
1606
+ @test sol. t[end ] == tspan[end ]
1607
+ @test sum (abs, sol. u[end ]) < 1
1608
+
1609
+ prob = ODEProblem {false} (lowered_dae_sys; u0_constructor = x -> SVector (x... ))
1610
+ @test prob. u0 isa SVector
1611
+ end
0 commit comments