Skip to content

Commit 50a4fc3

Browse files
author
fchen121
committed
Fix Linear transformation to diagonal system test
1 parent 22ab624 commit 50a4fc3

File tree

1 file changed

+30
-23
lines changed

1 file changed

+30
-23
lines changed

test/changeofvariables.jl

Lines changed: 30 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -62,36 +62,43 @@ new_sol = solve(new_prob, Tsit5())
6262
@test isapprox(sol[x][end], new_sol[x][end], rtol=1e-4)
6363

6464

65-
# # Linear transformation to diagonal system
66-
# @variables x(t)[1:3]
67-
# D = Differential(t)
68-
# A = [0. -1. 0.; -0.5 0.5 0.; 0. 0. -1.]
69-
# right = A.*transpose(x)
70-
# eqs = [D(x[1]) ~ sum(right[1, 1:3]), D(x[2]) ~ sum(right[2, 1:3]), D(x[3]) ~ sum(right[3, 1:3])]
65+
# Linear transformation to diagonal system
66+
@independent_variables t
67+
@variables x(t)[1:3]
68+
x = reshape(x, 3, 1)
69+
D = Differential(t)
70+
A = [0. -1. 0.; -0.5 0.5 0.; 0. 0. -1.]
71+
right = A*x
72+
eqs = vec(D.(x) .~ right)
7173

72-
# tspan = (0., 10.)
73-
# u0 = [x[1] => 1.0, x[2] => 2.0, x[3] => -1.0]
74+
tspan = (0., 10.)
75+
u0 = [x[1] => 1.0, x[2] => 2.0, x[3] => -1.0]
7476

75-
# @mtkcompile sys = System(eqs, t; defaults=u0)
76-
# prob = ODEProblem(sys,[],tspan)
77-
# sol = solve(prob, Tsit5())
77+
@mtkcompile sys = System(eqs, t; defaults=u0)
78+
prob = ODEProblem(sys,[],tspan)
79+
sol = solve(prob, Tsit5())
7880

79-
# T = eigen(A).vectors
81+
T = eigen(A).vectors
82+
T_inv = inv(T)
8083

81-
# @variables z(t)[1:3]
82-
# forward_subs = T \ x .=> z
83-
# backward_subs = x .=> T*z
84+
@variables z(t)[1:3]
85+
z = reshape(z, 3, 1)
86+
forward_subs = vec(T_inv*x .=> z)
87+
backward_subs = vec(x .=> T*z)
8488

85-
# new_sys = changeofvariables(sys, t, forward_subs, backward_subs; simplify=true)
89+
new_sys = changeofvariables(sys, t, forward_subs, backward_subs; simplify=true)
8690

87-
# new_prob = ODEProblem(new_sys, [], tspan)
88-
# new_sol = solve(new_prob, Tsit5())
91+
new_prob = ODEProblem(new_sys, [], tspan)
92+
new_sol = solve(new_prob, Tsit5())
8993

90-
# # test RHS
91-
# new_rhs = [eq.rhs for eq in equations(new_sys)]
92-
# new_A = Symbolics.value.(Symbolics.jacobian(new_rhs, z))
93-
# @test isapprox(diagm(eigen(A).values), new_A, rtol = 1e-10)
94-
# @test isapprox( new_sol[x[1],end], sol[x[1],end], rtol=1e-4)
94+
# test RHS
95+
new_rhs = [eq.rhs for eq in equations(new_sys)]
96+
new_A = Symbolics.value.(Symbolics.jacobian(new_rhs, z))
97+
A = diagm(eigen(A).values)
98+
A = sortslices(A, dims=1)
99+
new_A = sortslices(new_A, dims=1)
100+
@test isapprox(A, new_A, rtol = 1e-10)
101+
@test isapprox( new_sol[x[1],end], sol[x[1],end], rtol=1e-4)
95102

96103
# Change of variables for sde
97104
@independent_variables t

0 commit comments

Comments
 (0)