Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

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]




  • Prev by Date: Re: Preventing In-line Math Typesetting From Being Scaled Down in Text
  • Next by Date: Roots of a Jacobi polynomial
  • Previous by thread: Re: NDSolve in matrix form
  • Next by thread: default font