matrix - Julia 矩阵形式的 ODE 可能吗？

Question

我是 Julia 有前途的语言的新手，希望它能加速我的僵硬常微分方程。事情是这样的:1) 方程必须以矩阵形式定义，通过使用尺寸为 400x400 的质量、阻尼、刚度矩阵 outa Matlab。实现了二阶 ODE 的公共(public)状态空间表示。2)除了线性动力学之外，还有非线性力作用，这取决于它的某些状态。这些力必须在 ode 函数内定义。

然而，由于初始条件，状态变量根本不会改变，尽管应该改变。这是一个示例代码，具有较小的矩阵，用于原型(prototype)制作:

#Load packages
using LinearAlgebra
using OrdinaryDiffEq
using DifferentialEquations
using Plots

# Define constant matrices (here generated for the example)
const M=Matrix{Float64}(I, 4, 4) #  Mass matrix
const C=zeros(Float64, 4, 4) #  Damping matrix
const K=[10.0 0.0 0.0 0.0;  0.0 7.0 0.0 0.0; 0.0 0.0 6.0 0.0;0.0 0.0 5.0 0.0] #  Stiffness matrix

x0 = [0.0;0.0;0.0;0.0; 1.0; 1.0; 1.0; 1.0] #  Initial conditions
tspan = (0.,1.0) #  Simulation time span

#Define the underlying equation
function FourDOFoscillator(xdot,x,p,t)
  xdot=[-inv(M)*C -inv(M)*K; Matrix{Float64}(I, 4, 4) zeros(Float64, 4, 4)]*x
end

#Pass to Solvers
prob = ODEProblem(FourDOFoscillator,x0,tspan)
sol = solve(prob,alg_hints=[:stiff],reltol=1e-8,abstol=1e-8)
plot(sol)

我错过了什么？谢谢

参宿四

Answer 1

您没有改变输出，而是创建了一个新数组。如果您执行 xdot.=，它会起作用。

#Load packages
using LinearAlgebra
using OrdinaryDiffEq
using DifferentialEquations
using Plots

# Define constant matrices (here generated for the example)
const M=Matrix{Float64}(I, 4, 4) #  Mass matrix
const C=zeros(Float64, 4, 4) #  Damping matrix
const K=[10.0 0.0 0.0 0.0;  0.0 7.0 0.0 0.0; 0.0 0.0 6.0 0.0;0.0 0.0 5.0 0.0] #  Stiffness matrix

x0 = [0.0;0.0;0.0;0.0; 1.0; 1.0; 1.0; 1.0] #  Initial conditions
tspan = (0.,1.0) #  Simulation time span

#Define the underlying equation
function FourDOFoscillator(xdot,x,p,t)
  xdot.=[-inv(M)*C -inv(M)*K; Matrix{Float64}(I, 4, 4) zeros(Float64, 4, 4)]*x
end

#Pass to Solvers
prob = ODEProblem(FourDOFoscillator,x0,tspan)
sol = solve(prob,alg_hints=[:stiff],reltol=1e-8,abstol=1e-8)
plot(sol)