Re: NDSolve in matrix form
- To: mathgroup at smc.vnet.net
- Subject: [mg120491] Re: NDSolve in matrix form
- From: Glenn Carlson <g.crlsn at gmail.com>
- Date: Tue, 26 Jul 2011 07:05:35 -0400 (EDT)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
Here's a more complicated example.
The system is of the form A[i]'[t] + a[[i]] A[i][t] = Sum[G[t][[i]][[j]] A[j][t], {j,1,n}], where A is an array of solution functions with i = 1 to n. G[t][[]][[]] and a[[]] are specified arrays.
Note how I can work with different numbers of solution functions just by assigning a different value to n and defining the inputs G, a, AInit, and tf.
I'm still trying to make sense out why the single and double brackets have to be specified as they are, and why sometimes the index is after t and sometimes before t.
Thanks to all.
Glenn
-----
ClearAll;
n := 4
AInit := {1, 0, 0, 0};
G[t_] = {{0, t, 1, 0}, {0, 0, t, 1}, {1, 0, 0, t}, {t, 1, 0, 0}};
a = {5, 3, 2, 1};
tf := 5
(* Create a list of the solution functions excl. the argument t. *)
Asol = Table[A[i], {i, 1, n}];
(* Create a list of ODEs and ICs for NDSolve. *)
Table[Asol[[i]][t], {i, 1, n}];
eqns = Table[
Asol[[i]]'[t] + a[[i]] Asol[[i]][t] ==
Sum[G[t][[i]][[j]] Asol[[j]][t], {j, 1, n}], {i, 1, n}];
ic := Table[Asol[[i]][0] == AInit[[i]], {i, 1, n}];
sys = Join[eqns, ic];
Table[Asol[[i]]'[t] + a[[i]] Asol[[i]][t] ==
Sum[G[t][[i]][[j]] Asol[[j]][t], {j, 1, n}], {i, 1,
n}] // MatrixForm
Table[Asol[[i]][0] == AInit[[i]], {i, 1, n}] // MatrixForm
solution1 = NDSolve[sys, Asol, {t, 0, tf}, MaxSteps -> Infinity]
(* Create a list of solution functions incl. the argument t for ParametricPlot3D. *)
(*
fsolt=Table[A[i][t],{i,1,n}]
ParametricPlot3D[Evaluate[fsolt/.solution1],{t,0,tf},PlotPoints->1000,ColorFunction->(ColorData["Rainbow"][#4]&)]
*)
(* Plot solutions as functions of time. *)
(*
Plot[Evaluate[A[1][t]/.solution1],{t,0,tf},PlotRange->{0,2}]
Plot[Evaluate[A[2][t]/.solution1],{t,0,tf},PlotRange->{0,2}]
Plot[Evaluate[A[3][t]/.solution1],{t,0,tf},PlotRange->{0,2}]
Plot[Evaluate[A[4][t]/.solution1],{t,0,tf},PlotRange->{0,2}]
*)
Asolt = Table[A[i][t], {i, 1, n}];
Asum[t_] := Sum[(A[i][t] /. solution1), {i, 1, n}]
Plot[{Evaluate[Asolt /. solution1], Asum[t]}, {t, 0, tf}, PlotRange -> Automatic]