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, Jul 19 '17 at 15:07
For your math editing also explore the `amslatex` package, where you can find the `bmatrix` environment. — Lutz Lehmann, Jul 20 '17 at 14:09

score 5 · Accepted Answer · answered Jul 19 '17 at 14:55

5

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.

answered Jul 19 '17 at 14:55

Chris Rackauckas

18,645
3
50
81

Great. Thank you so much. I am impressed by the unique features of DifferentialEquations package. – Thanh Nguyen Jul 19 '17 at 16:07
Sorry, could you please give another example with matrices M, C, K (not scalars as used in your previous example). – Thanh Nguyen Jul 19 '17 at 16:47
It's the same if you make `M`, `C`, and `K` matrices. All that needs to be true is you have to update the values of `du`. Any mathematical statement that does that is valid. – Chris Rackauckas Jul 19 '17 at 17:03
I tried the following codes, but the results sol seems to be zeros for all values of t. Am I wrong somewhere here? `function eq(t,u,du); v=reshape([u...],Int(length(u)/2),2); du=reshape([v[:,2]; -M^(-1)*([10;10]-C*v[:,2]-K*v[:,1])],length(u)); end; M=[2.0 0.0; 0.0 1.0]; C=[0.0 0.0; 0.0 0.0]; K=[96.0 -32.0; -32.0 32.0]; u0=[0.0 ;0.0; 0.0; 0.0]; tspan=(0.0, 10.0); prob=ODEProblem(eq,u0,tspan); sol=solve(prob);` – Thanh Nguyen Jul 19 '17 at 17:18
You're not updating `du`, you're creating a new `du`. Notice that my code used `.=` instead of `=`. You should write it so you're updating `du` instead of replacing it. Or you can use the function signature `eq(t,u)` and return the array for the update (the MATLAB style), though this will take a hit to performance. – Chris Rackauckas Jul 19 '17 at 17:24
Even when I correct my code to update `du`, I still got zero `sol` for all `t`. Am I wrong somewhere? `function eq(t,u,du); v =reshape([u...],Int(length(u)/2),2); du .=reshape([v[:,2]; M^(-1)*([10;10]-C*v[:,2]-K*v[:,1])],length(u)); end; M=[2.0 0.0; 0.0 1.0]; C=[0.0 0.0; 0.0 0.0]; K=[96.0 -32.0; -32.0 32.0]; u0=[0.0 ;0.0; 0.0; 0.0]; tspan=(0.0, 10.0); prob=ODEProblem(eq,u0,tspan); sol=solve(prob);` – Thanh Nguyen Jul 19 '17 at 17:33
And here is my another attempt. This times, I got an error. `function eq(t,u,du) du .=[u[:,2]; M^(-1)*([10;10]-C*u[:,2]-K*u[:,1])]; end; M=[2.0 0.0; 0.0 1.0]; C=[0.0 0.0; 0.0 0.0]; K=[96.0 -32.0; -32.0 32.0]; u0=[0.0 0.0; 0.0 0.0]; tspan=(0.0, 10.0); prob=ODEProblem(eq,u0,tspan); sol=solve(prob);` – Thanh Nguyen Jul 19 '17 at 17:35
The matrix you're making doesn't look like it makes sense: the first row has a column vector, and the second row has a column vector? You likely meant to write that differently. Come to the chatroom to sort this out: https://gitter.im/JuliaDiffEq/Lobby – Chris Rackauckas Jul 19 '17 at 17:59

Using DifferentialEquations package in Julia to solve matrix ODE

1 Answers1