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 = 
  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.

  • 3
    $\begingroup$ Is there a reason why you have to multiply with the inverse on the right hand side rather than just multiplying with the original matrix on the left hand side? $\endgroup$ Commented Jan 17, 2018 at 14:59
  • $\begingroup$ I tried @SjoerdSmit's idea to put the original matrix on the left hand side of NDSolve. It seems to work fine. $\endgroup$
    – Chris K
    Commented Jan 17, 2018 at 18:38
  • $\begingroup$ I tried your second example. It gave a warning NDSolveValue::ivres, but otherwise produced a solution. Could you pinpoint the problem with this code? $\endgroup$
    – Chris K
    Commented Jan 18, 2018 at 21:28
  • $\begingroup$ Chris K, it will run, but the invmat, mat, val, etc. should change on every iteration of the loop, but they don't (as far as I know). The obvious way to do this is to make it a function, but if I do, the inverse evaluates symbolically, which takes forever for a 20x20. What I'm really trying to do is figure out how to tell mathematica that invmat MUST be populated with numbers before the inverse is evaluated on each iteration of NDSolve. $\endgroup$ Commented Jan 19, 2018 at 16:33

3 Answers 3


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}},
    {t, 0, 5}

Here are a couple plots:

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

enter image description here

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

enter image description here

  • $\begingroup$ Ah, very good; it does work. Mathematica complained about a list element that cannot be accessed, so I refrained from this. Using inv[x_List] := ... is the fix. $\endgroup$ Commented Jan 17, 2018 at 17:20
  • $\begingroup$ Thanks all for the feedback. In looking over the results, I think the example I provided was a bit too simple. I've added a more-detailed version of the code for clarity. $\endgroup$ Commented Jan 18, 2018 at 16:11

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.}]

enter image description here

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


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}]

Mathematica graphics

Looks like the other answers above.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.