# NDSolve problem with matrix valued function

I got a strange problem when using NDSolve to solve a matrix value function

σ[t] = {{σx[t], σxy[t]}, {σyx[t], σy[t]}}


when I try to solve it in a compact form.

Here is the initialization code:

Clear[t]
F = {{1, t}, {0, 1}};
LL = D[F, t].Inverse[F];
DD = 1/2 (LL + Transpose[LL]);
WW = 1/2 (LL - Transpose[LL]);
U = MatrixPower[Transpose[F].F, 1/2];
R = F.Inverse[U];
Ω = Simplify[D[R, t].Transpose[R]];
λ = 1;
μ = 1;


Now Ω is a 2 x 2 matrix with variable t. When I directly specify the form of σ, and I can get the problem solved:

σ = ({{σx[t], σxy[t]},{σyx[t], σy[t]}});
sol = NDSolve[{
D[σ, t] == (λ*Tr[DD])*IdentityMatrix + 2 μ*DD + Ω.σ + σ.Transpose[Ω],
(σ /. t -> 0) == ({{0, 0}, {0, 0}})
},
{σx, σy, σxy, σyx},
{t, 0, 5}]


However, if I try to define {σx[t], σxy[t]},{σyx[t], σy[t]}} as a single σ[t], which is more elegant, there will be problem:

Clear[σ];
func[X_?ArrayQ] := (λ*Tr[DD])*IdentityMatrix + 2 μ*DD + Ω.X + X.Transpose[Ω];
sol = NDSolve[{
D[σ[t], t] == func[σ][t]],
σ == ({{0, 0},{0, 0}})
}, σ, {t, 0, 5}]


The problem lays in the fact that the variable t in matrix Ω is not treated as the same as the t inside the NDSolve. Can you help me fix this?

Let everything that depend on t actually take t as argument, func especially:

Clear[t, F, LL, WW, U, Ω, λ, μ, R, DD]
F[t_] = {{1, t}, {0, 1}};
LL[t_] = D[F[t], t].Inverse[F[t]];
DD[t_] = 1/2 (LL[t] + Transpose[LL[t]]);
WW[t_] = 1/2 (LL[t] - Transpose[LL[t]]);
U[t_] = MatrixPower[Transpose[F[t]].F[t], 1/2];
R[t_] = F[t].Inverse[U[t]];
Ω[t_] = Simplify[D[R[t], t].Transpose[R[t]]];
λ = 1;
μ = 1;

Clear[σ];
func[X_?ArrayQ, t_] :=
(λ*Tr[DD[t]])*IdentityMatrix + 2 μ*DD[t] + Ω[t].X + X.Transpose[Ω[t]];
sol = NDSolve[{σ'[t] == func[σ[t], t], σ == {{0, 0}, {0, 0}}}
, σ
, {t, 0, 5}]

• It works, great! I first expect that I will need to use some Hold-like-command to treat the variable t in different parts of the code. Your code solves it all, thanks! Feb 16, 2013 at 19:25
• I have an additional question: Since σxy and σyx are identical, can we just tell Mathematica this and thus accelerate the calculation? Feb 20, 2013 at 21:59
• Add: I mean still using the compact way of expressing the equation. Feb 20, 2013 at 22:09