How to use DSolve with vectors without decomposing into vector components?

Question

Given a basic, two-dimensional ballistic trajectory problem, I can solve the equations of motion using DSolve (or NDSolve) by decomposing the vector equations of motion into equations with the scalar components, as follows:

DSolve[
  Flatten@{
    Thread[{rx''[t], ry''[t]} == {rx[t], -9.8 ry[t]}],
    Thread[{vx[t], vy[t]} == {rx'[t], ry'[t]}],
    {Thread[{rx[0], ry[0]} == {0, 0}],
     Thread[{vx[0], vy[0]} == {40 Cos[30 °], 
        40 Sin[30 °]}]}
    },
  {rx[t], ry[t], vx[t], vy[t]},
  t
  ];

But which general approaches are possible to do this just with "vectors"? (I'm not interested in having anyone solve this particular example per se.)

For example, something like this would be desirable:

DSolve[
 {
  Thread[m r''[t] == {0, -9.8} r[t]],
  Thread[v[t] == r'[t]],
  {Thread[r[0] == {0, 0}],
   Thread[v[0] == {40 Cos[30 °], 40 Sin[30 °]}]}
  },
 {r[t], v[t]},
 t
 ]

Surely this must be possible in some way? I understand that "vectors" are just Lists, and that Mathematica doesn't know that I want r[t], r'[t] or v[t] to be "vectors". However, using {r[t][[1]],r[t][[2]]} hasn't worked for me yet, either. I would also prefer to implement the method myself "manually" rather than rely on packages which I would have to see as a black box.

Answer 1

The idea (as in my comment) is to replace a vector-valued dependent variable by coordinate functions K[i] for i = 1, 2, ...; then call DSolve; and finally reassemble the solution coordinates into a vector-valued solutions.

ClearAll[reassembleSols, vectorize];
reassembleSols[sols_, fns_, var_, dim_] := Function[sol0, MapThread[
     Function[{f, s},
      If[FreeQ[f, var],
       f -> Function @@ {var, s},
       f -> s]],
     {fns, Partition[sol0, dim]}]
    ] /@ (Table[K[i][var], {i, 1 + Length[fns]*dim}] /. sols);

SetAttributes[vectorize, HoldAll];
vectorize[dim_, DSolve[ode_, f_, var_, rest___]] :=
 With[{fns = Flatten[{f}]},
  Block @@ 
   With[{depvars = Internal`ProcessEquations`FindDependentVariables[ode, var]},
    Hold[
     depvars,
     Module[{rules, sols},
      Do[(#[\[FormalT]_] = Table[K[i][\[FormalT]], {i, 1 + j*dim - dim, j*dim}]) &@
        depvars[[j]], {j, Length[depvars]}];
      rules = 
       Table[depvars[[j]] -> Table[K[i], {i, 1 + j*dim - dim, j*dim}],
        {j, Length[depvars]}];
      sols = DSolve[ode, Flatten[fns /. rules], var, rest];
      Clear @@ depvars;
      If[FreeQ[sols, DSolve],
       reassembleSols[sols, fns, var, dim],
       sols]
     ]]]

OP's example.

vectorize[2,
 DSolve[{m r''[t] == {{0, 0}, {0, -10}}. r[t], 
   v[t] == r'[t], {r[0] == {0, 0}, 
    v[0] == {40 Cos[30 °], 40 Sin[30 °]}}}, {r[t], 
   v[t]}, t]
 ]
(*
  {{r[t] -> {20 Sqrt[3] t,
      -I Sqrt[10] *
        E^(-((I Sqrt[10] t)/Sqrt[m])) (-1 + E^((2 I Sqrt[10] t)/Sqrt[m])) Sqrt[m]}, 
    v[t] -> {20 Sqrt[3], 
      10 E^(-((I Sqrt[10] t)/Sqrt[m])) (1 + E^((2 I Sqrt[10] t)/Sqrt[m]))}}}
*)

Note: The variables ({r, v} in the OP's example) are defined to be functions of the form

r[t_] = {K[1][t], K[2][t]}

inside the Block. These definitions cause the ode to evaluate to a componentwise system of equations when DSolve is called. When defined as functions this way, the derivatives of r are automatically computed properly in the system passed to DSolve. If we were to apply rules, this would not happen. Afterwards, these definitions are cleared in order to reassemble the solutions for r[t] etc. so that r[t] stays r[t] and does not evaluate to a vector of K's.