Tensorial functions

Question

I need to solve a system of ODEs for a function that has a complicated tensorial structure. I believe that NDSolve can deal with this scenario since version 9. Yet, in some cases it cannot figure out the dimensionality of the dependent variable from the initial conditions. That is why I want to map a tensorial function onto a tensor of functions. I do as follows:

Vector (rank-1):

Through[Array[y, {2}][t]]
(* {y[1][t], y[2][t]} *)

Matrix (rank-2):

Through /@ Through[Array[y, {2, 2}][t]]
(* {{y[1, 1][t], y[1, 2][t]}, {y[2, 1][t], y[2, 2][t]}} *)

Tensor (rank-3, naïve approach does not work):

Through /@ Through /@ Through[Array[y, {2, 2, 2}][t]]

Q. What is the proper way of automatically obtaining a rank-$k$ tensorial function?

Why is it needed?

The main goal is to minimize syntactic sugar for solving tensorial ODEs.

1. I will start from an example that can be found in the manual for NDSolve:

a={{1,  1,  1,  1},{1,  2,  1,  2},{1,  1,  3,  1},{1,  2,  1,  4}};
b={1,1,1,1};
NDSolve[{y''[x]+a.y[x]==0,y'[0]==y[0]==b},y,{x,0,8}]

This is a standard approach, where MA automatically figures out that y[x] is a vector. Nice.

2. But we can do the same slightly differently:

z[x_]=Through[Array[Y,{4}][x]];
s1=NDSolve[Join[Thread[z''[x]+a.z[x]==0],Thread[z'[0]==z[0]==b]],Array[Y,{4}],{x,0,8}];

MA solves here for 4 scalar functions instead {Y[1][x],Y[2][x],Y[3][x],Y[4][x]}, which is sometimes preferable. We can plot them as follows

Plot[Evaluate[z[x]/.s1],{x,0,8}]

3. Now, define a matrix function (two equivalent ways are possible)

w[x_]=Table[Y[i,j][x],{i,4},{j,4}]
w[x_]=Through/@Through[Array[Y,{4,4}][x]]

Notice that adding a dimension makes the syntax less and less transparent, and the commands become lengthier. Nonetheless, we can solve also this problem

s2=NDSolve[Join[Thread/@Thread[w''[x]+a.w[x]==0],Thread/@Thread[w'[0]==w[0]==a]]//Flatten,
    Array[Y,{4,4}]//Flatten,{x,0,8}];
Plot[Evaluate[w[x]/.s2],{x,0,8}]

What is disturbing is that going from a vector to a matrix equation more and more syntactic sugar needs to be added, like Thread/@Thread and Flatten. And it will be even less transparent for tensorial functions.

Extended question. What can be done to reduce this burden?

Answer 1

You can achieve what you want in the first part of your question using a Map and the levels argument. I've packed this into a function so you can change the function y, dimension, tensor-rank, or argument list (currently just {t}) if needs be:

genTensor[v_, dim_, rank_, args_] :=
  Map[Apply[#, args] &, Array[v, ConstantArray[dim, rank]], {rank}] 

(* 2d vector *)
genTensor[y,2,1,{t}]
(* result: {y[1][t], y[2][t]} *)

(* 2d matrix *)
genTensor[y,2,2,{t}]
(* result: {{y[1, 1][t], y[1, 2][t]}, {y[2, 1][t], y[2, 2][t]}} *)

(* 2d 3-tensor *)
genTensor[y,2,3,{t}]
(* result: {{{y[1, 1, 1][t], y[1, 1, 2][t]}, {y[1, 2, 1][t], y[1, 2, 2][t]}},
   {{y[2, 1, 1][t], y[2, 1, 2][t]},  {y[2, 2, 1][t], y[2, 2, 2][t]}}} *)