From: kogeddes@daisy.uwaterloo.ca (Keith O. Geddes)
Subject: Re: dumb maple question -- differential equations
Organization: University of Waterloo
Date: Thu, 17 Dec 1992 19:59:31 GMT

In article <DRCLARK.92Dec17110649@daisy.uwaterloo.ca> drclark@daisy.uwaterloo.ca (David R. Clark) writes:
>In article <BAGCHI.92Dec16184006@quip.eecs.umich.edu> bagchi@eecs.umich.edu (Ranjan Drzzzzt!  Bagchi) writes:
>
>   From: bagchi@eecs.umich.edu (Ranjan Drzzzzt!  Bagchi)
>   Subject: dumb maple question -- differential equations
>
>   ...   (header and sample session deleted)
>
>   There are several things I'd like to do with this data which 
>   I haven't had much luck with:
>
>   1)  I'd like to break v1(t), v2(t), pr(t) into separate functions.  While
>   this sounds trivial, the only way I've been able to do it is with 
>
>   V1 proc(t) op(2,op(2,pv(t))) end;
>
>   which looks fairly inelegant.
>
>Unfortunately, because Maple's sets are unordered, the hardcoded "2" in
>the first op call above will cause problems.  A more robust, but uglier,
>version is:
>V1 := proc(p)
>        local i,v;
>           v := pv(p);
>           for i in v do
>              if op(1,i) = 'v1(t)' then RETURN(op(2,i)) fi
>           od;
>           ERROR(`Could not find v1`);
>      end;


More to be recommended is the "trick" of substituting the solution
(a set of equations) into the desired variable name, as in:

V1 := T -> subs(pv(T), v1(t));

See the sample session below.

==============================================================================

    |\^/|     Maple V Release 2
._|\|   |/|_. Copyright (c) 1981-1992 by the University of Waterloo.
 \  MAPLE  /  All rights reserved. MAPLE is a registered trademark of
 <____ ____>  Waterloo Maple Software.
      |       Type ? for help.

> sys := { diff(v1(t),t) = a1*v1(t) - b1*v1(t)*pr(t),
>          diff(v2(t),t) = a2*v2(t) - b2*v2(t)*pr(t),
>           diff(pr(t),t) = -c*pr(t) + d * pr(t) * v1(t) + e *pr(t)*v2(t) };
                    d
          sys := {---- v1(t) = a1 v1(t) - b1 v1(t) pr(t),
                   dt

                d
              ---- v2(t) = a2 v2(t) - b2 v2(t) pr(t),
               dt

                d
              ---- pr(t) = - c pr(t) + d pr(t) v1(t) + e pr(t) v2(t)}
               dt

> params := {a1 = 1, b1 = .3,a2 = 2, b2 = .6, c = 3, d = .7, e = .5};
      params := {a1 = 1, b1 = .3, a2 = 2, b2 = .6, c = 3, d = .7, e = .5}

> sysinit := { v1(0) = 5, v2(0) = 5, pr(0) = 4 };
                  sysinit := {v1(0) = 5, v2(0) = 5, pr(0) = 4}

> eqns := subs(params, sys union sysinit);
                                            d
eqns := {v1(0) = 5, v2(0) = 5, pr(0) = 4, ---- v1(t) = v1(t) - .3 v1(t) pr(t),
                                           dt

      d
    ---- v2(t) = 2 v2(t) - .6 v2(t) pr(t),
     dt

      d
    ---- pr(t) = - 3 pr(t) + .7 v1(t) pr(t) + .5 v2(t) pr(t)}
     dt

> pv := dsolve(eqns, {v1(t),v2(t),pr(t)}, numeric):

> pv(1);
    {t = 1., v1(t) = 1.599162730, pr(t) = 5.059208376, v2(t) = .5114643376}

#
# 1)  I'd like to break v1(t), v2(t), pr(t) into separate functions.  While
# this sounds trivial, the only way I've been able to do it is with
#
# V1 :=  proc(t) op(2,op(2,pv(t))) end;
#
# which looks fairly inelegant.
#
#-----------------------------
# From: drclark@daisy.uwaterloo.ca (David R. Clark)
#
# Unfortunately, because Maple's sets are unordered, the hardcoded "2" in
# the first op call above will cause problems.  A more robust, but uglier,
# version is:
# V1 := proc(p)
#        local i,v;
#           v := pv(p);
#           for i in v do
#              if op(1,i) = 'v1(t)' then RETURN(op(2,i)) fi
#           od;
#           ERROR(`Could not find v1`);
#      end;
#-----------------------------
#
# More to be recommended is the "trick" of substituting the solution
# (a set of equations) into the desired variable name.
#
> V1 := T -> subs(pv(T), v1(t));
                         V1 := T -> subs(pv(T), v1(t))

> V1(1);
                                  1.599162730

> V1(1.5);
                                  1.579211333

> quit
Newsgroups: sci.math.symbolic