0
{a, b} = {6, 7};
t = 0;
dt = 0.001;
Do[{l, m} = With[{r = NSolve[{x - y == a, 2 x + y == b}, {x, y}]},
    Extract[r, Position[r, _?NumericQ]]];
    t = t + dt;
    a = a + (l + m)a dt;
    b = b + l m b dt, {1000}]

How to plot a(b), a(t) ?

Sektor
  • 3,320
  • 7
  • 27
  • 36
adam
  • 119
  • 1
  • You don't have functions a[t], b[t] produced by your code, just two individual numbers, a and b with values 2.45035 and 0.0609429, respectively. – murray Aug 20 '15 at 15:43

2 Answers2

7

Adam, I think this is probably what you mean:

{a, b} = {6, 7};
t = 0;
dt = 0.001;

loopvals = Part[#, 2, 1] &@
   Reap@
    Do[
     {l, m} = 
      With[{r = NSolve[{x - y == a, 2 x + y == b}, {x, y}]}, 
       Extract[r, Position[r, _?NumericQ]]];
     t = t + dt;
     a = a + (l + m) a dt;
     b = b + l m b dt;
     Sow[{a, b}],
     {1000}
    ];

results = Transpose@Prepend[loopvals, {6, 7}];
ListLinePlot[results, PlotLegends -> {"a", "b"}]

Mathematica graphics

Do loops return Null, so you have to actively save your values as they are calculated. Sow and Reap are one of the most efficient ways to do so, as you can see in this excellent tutorial from WRI on writing fast Mathematica code, and in countless posts on this site.

Your code has other problems as well.

  1. The system of equations you solve numerically at each iteration admits a simple analytical solution. You should really solve that symbolically once and then use the solution repeatedly.

For instance, if you solve the linear system, you obtain the following expressions for your recursion:

Clear[a, b, l, m]
{l, m} = {x, y} /. First@Solve[{x - y == a, 2 x + y == b}, {x, y}];
FullSimplify[{a + (l + m) a dt, b + l m b dt}]

(* Out: 
{a - 1/3 a (a - 2 b) dt, b + 1/9 b (-2 a + b) (a + b) dt}
*)
  1. Loops are typically a bad idea in Mathematica and they can be almost always avoided. You will want to look into alternative constructs. Read through this Q&A on this site: Alternatives to procedural loops and iterating over lists in Mathematica.

Combining both points, you can use the explicit recursion rule derived above and the idiomatic NestList command to apply this recursion rule repeatedly:

results = Transpose@
    NestList[
      Function[list,
        With[{a = list[[1]], b = list[[2]], dt = 0.001}, 
          {a - 1/3 a (a - 2 b) dt, b + 1/9 b (-2 a + b) (a + b) dt}]
      ],
      {6, 7},
      10000
    ];

Your original approach with NSolve took 8.34s to calculate 10,000 iterations; an approach still using Do, but using the pre-calculated solution (code not shown) takes 0.059s, yielding a 140x speedup; the idiomatic NestList approach above takes 0.002s to do the same, i.e. it is more than 4000x faster than the original, and 30x faster than Do.

Having said that, I do have one further comment: you need to learn to write better questions! Your question was barely intelligible, you had no code formatting, and you showed no effort to solve your problem yourself, and this is not the first question you asked on this site!

Community
  • 1
MarcoB
  • 67,153
  • 18
  • 91
  • 189
  • 2
    you are strict but just. someday, i'll ask a question which is as good as your answer, thank you. – adam Aug 20 '15 at 16:49
  • ...and if one wants to avoid component extractions: With[{dt = 0.001}, Transpose @ NestList[Function[list, list (1 + dt Total[list] {{-1/3, 2/3}, {-2/9, 1/9}}.list)], {6, 7}, 10]] – J. M.'s missing motivation Aug 20 '15 at 20:09
3

Just another way to do this:

m = {{1, -1}, {2, 1}};
fun = {{1 + (#1 + #2) 0.001, 0}, {0, 1 + #1 #2 0.001}} &;
nf[x_] := fun @@ (Inverse[m].x).x

So,

ListLinePlot[Transpose[NestList[nf, {6, 7}, 1000]],, PlotLegends -> {"a", "b"}]

yielding:

enter image description here

* Verbatim instructive comment from @Guesswhoitis.*

It's alright here, but in general, it's better to use LinearSolve[] instead of Inverse[]:

lsf = LinearSolve[{{1, -1}, {2, 1}}];
nf[x_] := fun @@ x.lsf[x].
ubpdqn
  • 60,617
  • 3
  • 59
  • 148