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

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

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;
(* 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]


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 8:04
Sorry, by k I mean c! –  Susan Aug 30 '13 at 9:07
@PlatoManiac +1 for visualizing the solution matrices ! Really nice :) ! –  Sektor Aug 30 '13 at 10:16

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]


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]


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