问题描述
我尝试在Julia中重现该示例,该示例在图中显示,取自Matrix Exponentiation
我向您展示我已经成功复制了Julia的练习。但我不知道如何针对感兴趣的范围引入向量t,例如t = -3:0.25:3.在矩阵中:[exp(u1 * t 0; 0 exp(u2 * t], u1 u2个特征值。
Julia>A=[0 1;1 0]
2×2 Array{Int64,2}:
0 1
1 0
F=eigen(A)
Eigen{Float64,Float64,Array{Float64,2},1}}
values:
2-element Array{Float64,1}:
-1.0
1.0
vectors:
2×2 Array{Float64,2}:
-0.707107 0.707107
0.707107 0.707107
D = diagm(exp.(F.values))
2×2 Array{Float64,2}:
0.367879 0.0
0.0 2.71828
P = F.vectors
13:06:08->>2×2 Array{Float64,2}:
-0.707107 0.707107
0.707107 0.707107
解决方法
在网站Fabian Dablander上 代码显示在R中实现该解决方案的代码。 这些是带给Julia的脚本:
using Plots
using LinearAlgebra
#Solving differential equations using matrix exponentials
A=[-0.20 -1;1 0] #[-0.40 -1;1 0.45] A=[0 1;1 0]
x0=[1 1]# [1 1] x0=[0.25 0.25] x0=[1 0]
tmax=20
n=1000
ts=LinRange(0,tmax,n)
x = Array{Float64}(undef,0)
x=x0
for i in 1:n
x=vcat(x,x0*exp(A*ts[i]))
end
plot(x)
plot(x[:,1],x[:,2])
#Solving differential equations finding eigenvalues and eigenvectors
A=[-0.20 -1;1 0] #A=[-0.40 -1;1 0.45] A=[0 1;1 0]
x0=[1 1]# [1 1] x0=[0.25 0.25] x0=[1 0]
tmax=20
n=500
# compute eigenvectors and eigenvalues
eig = eigen(A)
E =eig.vectors
λ=eig.values
# solve for the initial conditon
C =E\x0'
# create time steps
ts=LinRange(0,n,size(A,2))
for i in 1:n
x[i,:]=real.(C'*diagm(exp.(λ*ts[i]))*E)
end
plot(x)
plot(x[:,2])