How to pick out specific parts of the solution of a differential equation efficiently?

Question

I'm trying to solve a set of differential equations (in parallel) and then select specific parts of the solution. What I have at the moment (see below-I've simplified the differential equations, in fact they include matrices so that xx, yy and zz will be matrices) works, but uses a lot of memory. Ultimately, I only need the solutions of the differential equations at particular points- is there an efficient way to do this so that Mathematica only remembers the parts I need?

Thanks!

X = SparseArray[KroneckerProduct[{{0, 0, 0, 0, 0}, {0.5, 0, 0, 0, 0}, {1, 0, 0, 0,0}, {0, 0, 0, 0, 0}, {0, 1, 0.5, 0, 0}}, IdentityMatrix[3]]];
\[Rho]g = KroneckerProduct[{{0 , 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 0}, {0, 0, 0, 0, 1}}, {{1, 0, 0}, {0, 0, 0}, {0,
  0, 0}}];
sol = ParallelTable[NDSolve[{xx'[t] == xx[t], yy'[t] == c*yy[t], zz'[t] == zz[t],xx[0]==\[Rho]g, yy[0] == \[Rho]g, zz[0] == \[Rho]g}, {xx, yy, zz}, {t, 0, 2000}], {c, 0, 1000,10}];
a1[x_, i_] := (xx[t] /. sol[[i, 1]])
a2[x_, i_] := (yy[t] /. sol[[i, 1]])
a3[x_, i_] := (zz[t] /. sol[[i, 1]])
answer[t1_, i_] := Total[X.Transpose[a1[t1, i] - a2[t1, i] - a3[t1, i]], 2]
choosevalues = Table[answer[t1, i], {i, 1, 101}, {t1, 500 + (i - 1)*10, 500 + (i -1)*10 + 1000, 10}];

If I understand correctly it is a variation of this question Take a closer look at links provided in comments there. — Kuba, Aug 29 '13 at 10:00
Using something like this I'm not sure how to take into account the different values of c? Thanks! — Susan, Aug 29 '13 at 14:32

PlatoManiac · Answer 1 · 2013-08-30T11:56:55.147

Some food for brain! Will elaborate soon.

Clear[sol];
sol[tval_?NumericQ, cval_?NumericQ] := 
Module[{res, c}, 
  res = ParametricNDSolveValue[{xx'[t] == xx[t], yy'[t] == c*yy[t], 
  zz'[t] == zz[t], xx[0] == 2, yy[0] == 1, 
  zz[0] == 2}, {xx[tval] - yy[tval] - zz[tval]}, {t, 0, 100}, {c},
  MaxSteps -> 10^6, AccuracyGoal -> 2];
  res[cval]
  ];

This plays the role of your answer[t,i]

sol[15, 12]

{-1.53353*10^78}

Update:

First lets handle the equations properly. I inserted the HilbertMatrix to get nontrivial results. Your example matrix is too full of zeros.

var={xx[t],yy[t],zz[t]}=Table[Array[i[#1,#2][t]&,{15,15}],{i,{x,y,z}}];
{xx'[t],yy'[t],zz'[t]}=({xx[t],yy[t],zz[t]})/.a___[t]->a'[t]; 
X=SparseArray[KroneckerProduct[{{0,0,0,0,0},{0.5,0,0,0,0},{1,0,0,0,0},{0,0,0,0,0},  
 {0,1,0.5,0,0}},IdentityMatrix[3]]];
ρg=KroneckerProduct[HilbertMatrix[5],{{1,0,0},{0,0,0},{0,0,0}}];
A = RandomReal[{-1, 1}, {15, 15}];
B = c RandomReal[{-.5, .5}, {15, 15}]; (* B Depends on c *)
Com[x_, y_] := x.y - y.x;
eqs=Flatten@(Join@@MapThread[#1==#2&,#,2]&/@{
   (* Equations {LHS matrix, RHS matrix} *)
   (* the xx'[t]==Com[(A-B[c]),xx[t]] *)
   {xx'[t],Com[(A - B), xx[t]]},
   {yy'[t],c yy[t]},
   {zz'[t],zz[t]},
   (* Equations LHS\[Equal]RHS with {LHS matrix & RHS matrix} *)
   {xx[t],ρg}/.t-> 0,
   {yy[t],ρg}/.t-> 0,
   {zz[t],ρg}/.t-> 0});

Now the version 8 compatible function sol8 that takes different values of c. The process can be made more optimal. Look at the suggested post by @Kuba.

Clear[sol8];
sol8[tval_?NumericQ, cval_?NumericQ] := Module[{rule},
   rule =Dispatch@Flatten@NDSolve[
         Evaluate[eqs /. c -> cval], 
         Flatten@var /. a___[t] -> a,
         {t, 0, 20},AccuracyGoal -> 2];
  (X . (xx[t] - yy[t] - zz[t])) /. rule /.t -> tval];

Now testing!

DistributeDefinitions[sol8];
choosevalues = ParallelTable[Evaluate@sol8[tval, -Cos@cval], 
 {cval, 1, 10, 1}, {tval, 1, 20, .5}];

Check the dimension of choosevalues

Dimensions@choosevalues

{10, 39, 15, 15}

In Mathematica 9 you can use Image3D to visualize the solution matrices as voxels.

Image3D[choosevalues]

enter image description here

Hope this gets you going!

BR

OK, thanks! The only problem I might have is with using ParametricNDSolveValue, since I only have Mathematica 8. Do you know if there's an alternative? — Susan, Aug 29 '13 at 11:11
I'm trying to solve the differential equations where xx, yy and zz are all matrices (with dimensions 15x15). Sorry, perhaps I haven't been very clear- the differential equations are really more complicated and include multiplication by other matrices etc. — Susan, Aug 29 '13 at 11:31
Thanks for your answer! My differential equations are a bit more complicated, for example something like k*A.xx[t], where A is another 15x15 matrix. I've tried putting this into your answer and I get the error message 'Dot::dotsh: Tensors...and...have incompatible shapes' for the 'eqs=...' line. Please can you explain a bit more what this line is doing, so that I can see where it's going wrong? Thanks! — Susan, Aug 30 '13 at 08:04
@PlatoManiac +1 for visualizing the solution matrices ! Really nice :) ! — Sektor, Aug 30 '13 at 10:16
@PlatoManiac Thanks for this! I'm still a bit confused though as to why you write the equations as something like: (c (RandomReal[{-1,1},{15,15}]) . xx[t])? Why do you not need to specify that it needs to be multiplied by c? If instead, I wanted the equation to be something like Com[(A-B[c]),xx[t]] where I've already specified B (a 15x15 matrix, which depends on c) and the general function Com[x_, y_] := x.y - y.x; how would I write the equations? — Susan, Aug 30 '13 at 10:32

score 1 · Answer 2 · answered Aug 29 '13 at 10:23

Yes, there is an efficient way. It is important to see that the range of solution does not have to start from where initial conditions are. So you can request a solution to be between any points, far removed from the initial conditions. Here is an example:

ic = y[0] == 1;
ode = y'[x] == y[x] Cos[x + y[x]];
sol = NDSolve[{ode, ic}, y, {x, 0, 30}]
Plot[y[x] /. sol, {x, 0, 30}, PlotRange -> All, Evaluated -> True]

Mathematica graphics

But suppose you want the solution only from x=15 to x=16 ? You can just ask for that part, and still use the same initial conditions at x=0

ic = y[0] == 1;
ode = y'[x] == y[x] Cos[x + y[x]];
sol = NDSolve[{ode, ic}, y, {x, 15, 16}]
Plot[Evaluate[y[x] /. sol], {x, 15, 16},PlotRange -> {{0, 30}, {0, 1}}, PlotStyle -> Red]

Mathematica graphics

You can even do more. You can as NDSolve to solve for small segments, and treat it as an integration step itself, and use the new solution to solve for the next segment.

Sorry, I've made a mistake with the initial conditions (which should be matrices), I've changed this in the question.The problem is that I want the solution over as wide a range of values as possible, but I only need specific values (at a given interval) if that makes sense? — Susan, Aug 29 '13 at 10:34

How to pick out specific parts of the solution of a differential equation efficiently?

2 Answers2