Dynamically updating parameter each iteration of NDSolve

Question

I am struggling to update the value of a parameter inside NDSolve that changes with the current value of my state. See an extremely simplified version of the code below.

x[t] = {x1[t], x2[t]};

x'[t] = D[x[t], t];

Cmat = {{10, x1[t]}, {x2[t], 30}};

Cmatinv = Inverse[Cmat /. x1[t] -> 3 /. x2[t] -> 5]

{{2/19, -(1/95)}, {-(1/57), 2/57}}

NDSolve[{x'[t] == Cmatinv. x[t], x1[0] == 3, x2[0] == 5}, 
   x[t], {t, 0, 5}];

Essentially what I am trying to do is recalculate Cmatinv at each iteration within NDSolve. The problem is that, in my real code, Cmat is a 20x20 with large expressions for each element, so I can't do a symbolic inverse like I could for the simplified code above. To initialize the differential equation, I substitute the initial conditions into Cmat so that the inverse is done on a numeric matrix and is therefore extremely fast. I can't figure out how to then update this at each iteration with the new states within NDSolve. The procedure I want to accomplish within each NDSolve iteration is: calculate x[t], calculate Cmat (numeric), calculate Cmatinv, solve differential equation. I've been troubleshooting this problem for months, so any help would be appreciated. I've seen similar questions answered, but those evaluate the relevant parameters as functions of the state. As far as I can tell, I can't do this, because the inverse would end up being symbolic.This is my first time posting, so I apologize if I'm missing any rules/information.

I am including a more detailed version of my code, because the simplified version doesn't really show exactly what I'm struggling with.

(*Physical Parameters*)
param1 = 8.99*10^9;
param2 = 9.81;

(*System Parameters*)
n = 20;
tmax = 0.1;

(*Spring Coefficient*)
springparam = 2.9*10^-6;

(*Material Properties and Mass Calculation*)
length = 0.06527392544827149;
l = length/n;
m = 8.668377299530452*10^-8;

(*Field Calculations*)
sysparam = 6150;
sphparam = ConstantArray[sysparam, n];
fieldparam = {60531.496062992126, 0, 0};

(*System Setup*)
h[1] = {0, 0, 0};
h[i_] := h[i - 1] + 
  l {Sin[Subscript[θ, i - 1][t]], -Cos[Subscript[θ, i - 1][t]], 0}
hinge = Table[h[ii], {ii, n}];
com = Table[
     hinge[[i]] + l/2 {Sin[Subscript[θ, i][t]], -Cos[Subscript[θ, i][t]], 0}
   , {i, n}];
vcom = Table[D[com[[i]], t], {i, n}];
comypos = Table[-(2 i - 1) l/2, {i, n}];

initpos = Table[Subscript[θ, i][t] -> 0, {i, n}];

r[i_, j_] := Sqrt[(com[[i, 1]] - com[[j, 1]])^2 + (com[[i, 2]] - com[[j, 2]])^2]
r[i_, i_] := 0.0009242787843475243;
invmat = (Table[(r[i, j])^-1, {i, n}, {j, n}]);

mat = 1/param1 Inverse[invmat /. initpos];
val = mat.sphparam;

(*Generalized Force Setup*)
forceparam = Table[fieldparam*val[[i]], {i, n}];
genforceparam = Table[Sum[
    forceparam[[k]].D[com[[k]], Subscript[θ, i][t]]
  , {k, n}], {i, n}];

(*Langrangian Setup*)
u = Table[
   m param2 (com[[i, 2]] - comypos[[i]]) + 
    1/2 springparam*(Subscript[θ, i][t] - Subscript[θ, i - 1][t])^2, {i, 1, n}];
u[[1]] = m param2 (com[[1, 2]] - comypos[[1]])
  + 1/2 springparam (Subscript[θ, 1][t] - π/2)^2;
k = Table[1/2 m vcom[[i]].vcom[[i]] + 1/24 m l^2 Subscript[θ, i]'[t]^2, {i, n}];

L = Sum [k[[i]], {i, n}] - Sum[u[[i]], {i, n}];
eulerlagrange[i_] := 
  D[D[L, Subscript[θ, i]'[t]], t] - D[L, Subscript[θ, i][t]] == genforceparam[[i]];
diffeq = Join[Table[eulerlagrange[i], {i, n}], 
   Table[Subscript[θ, i][0] == 0, {i, n}], 
   Table[Subscript[θ, i]'[0] == 0, {i, n}]];

evalel = 
 NDSolveValue[diffeq, 
  Join[Table[Subscript[θ, i][t], {i, n}], 
   Table[Subscript[θ, i]'[t], {i, n}]], {t, 10^-6, tmax }, 
  Method -> {"EquationSimplification" -> "Residual"}]

You can see here that, in order to get my generalized force for the Euler-Lagrange equation, I have to solve for val, meaning that I need to invert invmat which depends on com.

Answer 1

Since NDSolveValue works with vectors, I would tweak @HenrikSchumacher's answer as follows:

inv[x_List] := LinearSolve[{{10, x[[1]]}, {x[[2]], 30}}, x]

sol = NDSolveValue[
    {x'[t] == inv[x[t]], x[0] == {3, 5}},
    x,
    {t, 0, 5}
];

Here are a couple plots:

Plot[sol[t], {t, 0, 5}]

ParametricPlot3D[{t, sol[t][[1]], sol[t][[2]]}, {t, 0, 5}]

Answer 2

How about this?

ClearAll[x, y];
F[t_, x_, y_] := LinearSolve[{{10., x}, {y, 30.}}, {x, y}];
{x, y} = NDSolveValue[
  {
   {x'[t], y'[t]} == F[t, x[t], y[t]], 
   x[0] == 3, 
   y[0] == 5
  }, 
  {x, y}, 
  {t, 0, 5}
]
ParametricPlot3D[{t, x[t], y[t]}, {t, 0., 5.}]

Side remark: Avoid Inverse for all matrices bigger than, say $6 \times 6$.

Answer 3

To follow-up on my comment above, on @SjoerdSmit's observation that you could multiply the left hand side by the original matrix to avoid inverting it:

sol = NDSolve[{{{10, x1[t]}, {x2[t], 30}}.{x1'[t], x2'[t]} == {x1[t], x2[t]},
  x1[0] == 3, x2[0] == 5}, {x1, x2}, {t, 0, 5}];
Plot[{x1[t], x2[t]} /. sol, {t, 0, 5}]

Looks like the other answers above.