Linear ODEs in Julia


# Package to integrate ODEs

using DifferentialEquations 

# Package for Linear Algebra. We will use it to 
# compute the identiy matrix and other quantities.

using LinearAlgebra

# Taken from DifferentialEquations documentation
# Take a 4x4 matrix
C= Matrix(rand(2,2))
A  =C - transpose(C) # A is antisymmetric
Φ = [1.0 0.0
    0.0 1.0]  # Initial condition

2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0

# We will now set up the ODE

# Define a time interval. It is not necessary  
# to indicate a set of values where solutions will be
# computed.
tspan = (0.0,10)
f(X,p,t) = A*X # Define the ODE as a simple matrix product.
prob = ODEProblem(f,Φ,tspan) # We now set the ODE problem
sol = solve(prob) # and use the interface for solutions

retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 9-element Vector{Float64}:
  0.0
  0.007524128284904676
  0.08276541113395143
  0.6014725571134181
  1.8105220007270149
  3.624653687906628
  5.929465071505989
  8.86517836071115
 10.0
u: 9-element Vector{Matrix{Float64}}:
 [1.0 0.0; 0.0 1.0]
 [0.9999995009985435 -0.0009990008328333417; 0.0009990008328333417 0.9999995009985435]
 [0.9999396214263468 -0.010988789821184267; 0.010988789821184267 0.9999396214263468]
 [0.9968129397132733 -0.07977445212055055; 0.07977445212055055 0.9968129397132733]
 [0.9712455730519747 -0.23807989311744; 0.23807989311744 0.9712455730519747]
 [0.886414169475009 -0.4628929750011257; 0.4628929750011257 0.886414169475009]
 [0.7057799028763244 -0.7084311352335161; 0.7084311352335161 0.7057799028763244]
 [0.3836446175188135 -0.9234807265626944; 0.9234807265626944 0.3836446175188135]
 [0.24067966299298466 -0.970604524919966; 0.970604524919966 0.24067966299298466]

# sol is a function which gives the interpolated solution

D=sol(pi/4)
transpose(D)*D #D is in SO(2,R)

2×2 Matrix{Float64}:
 1.0  0.0
 0.0  1.0


# ## Plotting the solutions

using Plots 

plot(sol)


# Let us check Liouville theorem

det_of_solution(t)= det(sol(t))
plot(det_of_solution,0,10,label="Determinant of Solution")


# #Idea of Geometric integration
Tmax=1e3
tspan = (0.0,Tmax)
prob = ODEProblem(f,Φ,tspan) 
sol = solve(prob)
plot(t->opnorm(sol(t)),0,Tmax, label="Matrix Norm of Solution")

2×2 Matrix{Float64}:
 0.236086  0.281035
 0.413808  0.859431


# ## References
# - [https://docs.sciml.ai/DiffEqDocs/stable/types/nonautonomous_linear_ode/]
# - [https://en.wikipedia.org/wiki/Geometric_integrator]

#N=4
#C=randn(N,N)
C = [-4. -17.; 2. 2.]; N=2
eigvals(C)

2-element Vector{ComplexF64}:
 -1.0 - 5.0im
 -1.0 + 5.0im

rmax=maximum(real(eigvals(C)))

-1.0

A=C-(1+rmax)*Matrix{Float64}(I, N, N)
eigvals(A)

tspan = (0.0,10)
f(X,p,t) = A*X # Define the ODE as a simple matrix product.
prob = ODEProblem(f,Matrix(I,N,N),tspan) # We now set the ODE problem
sol = solve(prob) # and use the interface for solutions

retcode: Success
Interpolation: specialized 4th order "free" interpolation, specialized 2nd order "free" stiffness-aware interpolation
t: 67-element Vector{Float64}:
  0.0
  8.253669896410068e-5
  0.0009079036886051074
  0.009161573585015174
  0.03642532985991913
  0.08322657363638775
  0.14405347840160138
  0.21873814971959693
  0.3007827769847833
  0.4011529239453252
  ⋮
  8.245679407524333
  8.439735737708515
  8.65362921925336
  8.869802813007224
  9.078169267000215
  9.310206259925518
  9.545976184087746
  9.773112479875305
 10.0
u: 67-element Vector{Matrix{Float64}}:
 [1.0 0.0; 0.0 1.0]
 [0.9996697919065034 -0.001403008038131507; 0.000165059769191942 1.0001649712140792]
 [0.9963609841109974 -0.015420303188974556; 0.0018141533163499473 1.0018034440600472]
 [0.9626162914807793 -0.1542724155589438; 0.01814969594811103 1.0170653793251123]
 [0.8434967778414431 -0.5937853808643483; 0.06985710363109979 1.0530680887347423]
 [0.6184491627947923 -1.2646161266984621; 0.14877836784687787 1.0647842663354254]
 [0.3081338874135226 -1.941728284206627; 0.2284386216713677 0.9934497524276255]
 [-0.05929001871454903 -2.4269186575837383; 0.28551984206867476 0.797269507491475]
 [-0.3936777303996087 -2.5111815231337737; 0.29543312036867914 0.4926216307064276]
 [-0.6464542928516797 -2.0644782303341715; 0.242879791804021 0.0821850825603812]
 ⋮
 [-0.00018225355195093502 0.00033838844568805386; -3.9810405375064815e-5 -0.00030168476807613067]
 [8.149191612559598e-5 0.0007164332456619066; -8.428626419551758e-5 -0.0001713668764609588]
 [0.0001998745686797082 0.0003859151942206503; -4.540178755536969e-5 6.366920601359865e-5]
 [0.0001001677530307758 -0.0001723103667108356; 2.02718078483339e-5 0.0001609831765757786]
 [-4.942515048368562e-5 -0.0003817167188956619; 4.4907849281842276e-5 8.529839736184292e-5]
 [-0.000105045563429498 -0.00016489405508461917; 1.9399300598189595e-5 -4.684766163492794e-5]
 [-3.375919759609567e-5 0.0001388557406014767; -1.6335969482526954e-5 -8.276710604367908e-5]
 [4.3444416309737225e-5 0.0001896802838342912; -2.231532750991602e-5 -2.3501566220013033e-5]
 [5.058194001241343e-5 3.929918938196095e-5; -4.623434044936121e-6 3.671163787760482e-5]

plot(t->opnorm(sol(t)),0,10, label="Matrix Norm of Solution")
plot!(t->exp(rmax*t),0,10, label="e^{-rt}")

SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
U factor:
2×2 Matrix{Float64}:
  0.990519  -0.137379
 -0.137379  -0.990519
singular values:
2-element Vector{Float64}:
 17.630233797973155
  1.4747393765696442
Vt factor:
2×2 Matrix{Float64}:
 -0.240316  -0.970695
 -0.970695   0.240316

one(3)