Using DifferentialEquations package in Julia to solve matrix ODE

Question

I would like to solve:

[\mathbf{M} \ddot{ \mathbf{U} }+ \mathbf{C} \dot{ \mathbf{U} }+ \mathbf{K} \mathbf{U} = \mathbf{P}(t)]

Or, in state-space form:

[\dot{\mathbf{Y}}=f(\mathbf{Y},t)]

where:

[\mathbf{Y} = \left[ \begin{array}{ c} \mathbf{U} \ \dot{ \mathbf{U} \end{array} \right] ]

and:

[f( \mathbf{Y} ,t)= \left[ \begin{array}{ c} \dot{ \mathbf{U} }\ \mathbf{M}^{-1} \mathbf{P} (t)- \mathbf{M} ^{-1} \mathbf{C} \dot{ \mathbf{U} }- \mathbf{M} ^{-1} \mathbf{K} \mathbf{U} \end{array} \right] ]

I tried the following code in Julia, using

\mathbf{M} = \left[ \begin{array}{ cc} 2&0\ 0&1 \end{array} \right];

\mathbf{C} = \left[ \begin{array}{ cc} 0&0\ 0& 0 \end{array} \right];

\mathbf{K} = \left[ \begin{array}{ cc} 96&-32\ -32& 32 \end{array} \right];

\mathbf{P} (t)= \left[ \begin{array}{ c} 10\ 10 \end{array} \right]

.

using DifferentialEquations
function eq(t,u,du)
    v=reshape([u...],Int(length(u)/2),2)
    du=reshape([v[:,2];-[10;10]-M\C*v[:,2]-M\K*v[:,1]],length(u))
end
u0=[0;0;0;0];
tspan=(0.0,10.0);
prob=ODEProblem(eq,u0,tspan)
sol=solve(prob)

But running these lines of code results to this error:

ERROR: InexactError()

I am using Julia ver. 0.5.2.

Please help me. Thank you.

Your explicit equation should read U''=M^(-1)*(P(t)-C*U'-K*U). In your version the M^(-1) before P(t) is missing, both in the transformed equations and in the code. -- You should have defined Y, not \dot Y as [U;U']. — Lutz Lehmann
For your math editing also explore the amslatex package, where you can find the bmatrix environment. — Lutz Lehmann

Chris Rackauckas Chris Rackauckas · Accepted Answer · 2017-07-19T14:55:53

Your problem is that DifferentialEquations.jl respects your input type.

u0=[0;0;0;0];

This is an array of integers, and so this means that your problem will be evolving an array of integers. In the first step, it finds out that it the calculation returns floating point numbers, and so it doesn't know what to make u because you said it had to be an array of integers. The fix for that is to say that your problem is on floating point numbers:

u0=[0.0;0.0;0.0;0.0];

Now that evolves correctly.

But let's go another step. DifferentialEquations.jl respects your input type, so DifferentialEquations.jl will make it into a problem on matrices just by making the initial condition a matrix. No reshaping is necessary if you make u0 a matrix:

u0=[0.0 0.0
    0.0 0.0]

Once you do that, you just write an equation which directly works on matrices. Example:

using DifferentialEquations
M = 1
K = 1
C = 1
function eq(t,u,du)
  du .= u[:,2] .- 10 - M./C.*u[:,2] - M.\K.*u[:,1]
end
u0=[0.0 0.0
    0.0 0.0]
tspan=(0.0,10.0);
prob=ODEProblem(eq,u0,tspan)
sol=solve(prob)

I'm not sure I got the equation you were trying to solve correct because it's extremely hard to read, but that should get you very close to what you want.

Using DifferentialEquations package in Julia to solve matrix ODE

1 Answers