ODE in Julia -matrix form possible?

Question

I am new to the promising language of Julia, in the hope that it will accelerate my stiff ordinary differential equations. Here is the thing: 1) The equation must be defined in matrix form, by using mass, damping, stiffness matrices outa Matlab with dimension 400x400. The common state-space represenation for 2nd order ODE's is implemented. 2) Apart from the linear dynamics, there are nonlinear forces acting, which depend on certain states of it. These forces must be defined inside the ode function.

However the state variables do not change at all, although the should, due to the inital conditions. Here is an example code, with smaller matrices, for prototyping:

#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)

What am I missing? Thanks

Betelgeuse

You're not mutating the output, and instead creating a new array. If you do xdot.= it should work. — Chris Rackauckas

Chris Rackauckas Chris Rackauckas · Accepted Answer · 2020-01-07T13:40:01

You're not mutating the output, and instead creating a new array. If you do xdot.= it works.

#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)

ODE in Julia -matrix form possible?

1 Answers