Блог

Главная - Вопросы по программированию - Реализовать решение системы дифференциальных уравнений с использованием экспоненциальной матрицы в Julia

Реализовать решение системы дифференциальных уравнений с использованием экспоненциальной матрицы в Julia

#arrays #matrix #julia #differential-equations #exponentiation

#массивы #матрица #julia #дифференциальные уравнения #возведение в степень

Вопрос:

Я пытаюсь воспроизвести пример в Julia, который я показываю на рисунке и взят из возведения матрицы в степень
введите описание изображения здесь

Я покажу вам, как далеко мне удалось воспроизвести упражнение в 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},Array{Float64,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

Комментарии:

1. Просто используйте exp(A*t)*u0 . Больше ничего не нужно.

2. Обратите внимание, что если вы посмотрите, что он делает, он в конечном итоге делает в основном то, что вы бы написали в любом случае

3. возможно ли улучшить код?

4. C=[0 1;1 0] x0=[1 0] ∆t = .25 T = exp(C*∆t) x = x0 x1 =x0 for i = 1:100 x = x*T #repeatedly multiply by T x1=vcat(x1, x) # amp; store current x1(t) in the array x1 end

Ответ №1:

На сайте Fabian Dablander code показаны коды на 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, 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 =Ex0'
# create time steps
ts=LinRange(0,tmax,n)
x = Array{Float64}(undef,n, size(A,2))
for i in 1:n
   x[i,:]=real.(C'*diagm(exp.(λ*ts[i]))*E)
 end
plot(x)
plot(x[:,1],x[:,2])