home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #30 / NN_1992_30.iso / spool / sci / math / symbolic / 3308 < prev    next >
Encoding:
Text File  |  1992-12-17  |  4.4 KB  |  134 lines
  1. Newsgroups: sci.math.symbolic
  2. Path: sparky!uunet!zaphod.mps.ohio-state.edu!cs.utexas.edu!torn!watserv2.uwaterloo.ca!watdragon.uwaterloo.ca!daisy.uwaterloo.ca!kogeddes
  3. From: kogeddes@daisy.uwaterloo.ca (Keith O. Geddes)
  4. Subject: Re: dumb maple question -- differential equations
  5. Message-ID: <BzF678.M3J@watdragon.uwaterloo.ca>
  6. Sender: news@watdragon.uwaterloo.ca (USENET News System)
  7. Organization: University of Waterloo
  8. References: <BAGCHI.92Dec16184006@quip.eecs.umich.edu> <DRCLARK.92Dec17110649@daisy.uwaterloo.ca>
  9. Date: Thu, 17 Dec 1992 19:59:31 GMT
  10. Lines: 122
  11.  
  12. In article <DRCLARK.92Dec17110649@daisy.uwaterloo.ca> drclark@daisy.uwaterloo.ca (David R. Clark) writes:
  13. >In article <BAGCHI.92Dec16184006@quip.eecs.umich.edu> bagchi@eecs.umich.edu (Ranjan Drzzzzt!  Bagchi) writes:
  14. >
  15. >   From: bagchi@eecs.umich.edu (Ranjan Drzzzzt!  Bagchi)
  16. >   Subject: dumb maple question -- differential equations
  17. >
  18. >   ...   (header and sample session deleted)
  19. >
  20. >   There are several things I'd like to do with this data which 
  21. >   I haven't had much luck with:
  22. >
  23. >   1)  I'd like to break v1(t), v2(t), pr(t) into separate functions.  While
  24. >   this sounds trivial, the only way I've been able to do it is with 
  25. >
  26. >   V1 proc(t) op(2,op(2,pv(t))) end;
  27. >
  28. >   which looks fairly inelegant.
  29. >
  30. >Unfortunately, because Maple's sets are unordered, the hardcoded "2" in
  31. >the first op call above will cause problems.  A more robust, but uglier,
  32. >version is:
  33. >V1 := proc(p)
  34. >        local i,v;
  35. >           v := pv(p);
  36. >           for i in v do
  37. >              if op(1,i) = 'v1(t)' then RETURN(op(2,i)) fi
  38. >           od;
  39. >           ERROR(`Could not find v1`);
  40. >      end;
  41.  
  42.  
  43. More to be recommended is the "trick" of substituting the solution
  44. (a set of equations) into the desired variable name, as in:
  45.  
  46. V1 := T -> subs(pv(T), v1(t));
  47.  
  48. See the sample session below.
  49.  
  50. ==============================================================================
  51.  
  52.     |\^/|     Maple V Release 2
  53. ._|\|   |/|_. Copyright (c) 1981-1992 by the University of Waterloo.
  54.  \  MAPLE  /  All rights reserved. MAPLE is a registered trademark of
  55.  <____ ____>  Waterloo Maple Software.
  56.       |       Type ? for help.
  57.  
  58. > sys := { diff(v1(t),t) = a1*v1(t) - b1*v1(t)*pr(t),
  59. >          diff(v2(t),t) = a2*v2(t) - b2*v2(t)*pr(t),
  60. >           diff(pr(t),t) = -c*pr(t) + d * pr(t) * v1(t) + e *pr(t)*v2(t) };
  61.                     d
  62.           sys := {---- v1(t) = a1 v1(t) - b1 v1(t) pr(t),
  63.                    dt
  64.  
  65.                 d
  66.               ---- v2(t) = a2 v2(t) - b2 v2(t) pr(t),
  67.                dt
  68.  
  69.                 d
  70.               ---- pr(t) = - c pr(t) + d pr(t) v1(t) + e pr(t) v2(t)}
  71.                dt
  72.  
  73. > params := {a1 = 1, b1 = .3,a2 = 2, b2 = .6, c = 3, d = .7, e = .5};
  74.       params := {a1 = 1, b1 = .3, a2 = 2, b2 = .6, c = 3, d = .7, e = .5}
  75.  
  76. > sysinit := { v1(0) = 5, v2(0) = 5, pr(0) = 4 };
  77.                   sysinit := {v1(0) = 5, v2(0) = 5, pr(0) = 4}
  78.  
  79. > eqns := subs(params, sys union sysinit);
  80.                                             d
  81. eqns := {v1(0) = 5, v2(0) = 5, pr(0) = 4, ---- v1(t) = v1(t) - .3 v1(t) pr(t),
  82.                                            dt
  83.  
  84.       d
  85.     ---- v2(t) = 2 v2(t) - .6 v2(t) pr(t),
  86.      dt
  87.  
  88.       d
  89.     ---- pr(t) = - 3 pr(t) + .7 v1(t) pr(t) + .5 v2(t) pr(t)}
  90.      dt
  91.  
  92. > pv := dsolve(eqns, {v1(t),v2(t),pr(t)}, numeric):
  93.  
  94. > pv(1);
  95.     {t = 1., v1(t) = 1.599162730, pr(t) = 5.059208376, v2(t) = .5114643376}
  96.  
  97. #
  98. # 1)  I'd like to break v1(t), v2(t), pr(t) into separate functions.  While
  99. # this sounds trivial, the only way I've been able to do it is with
  100. #
  101. # V1 :=  proc(t) op(2,op(2,pv(t))) end;
  102. #
  103. # which looks fairly inelegant.
  104. #
  105. #-----------------------------
  106. # From: drclark@daisy.uwaterloo.ca (David R. Clark)
  107. #
  108. # Unfortunately, because Maple's sets are unordered, the hardcoded "2" in
  109. # the first op call above will cause problems.  A more robust, but uglier,
  110. # version is:
  111. # V1 := proc(p)
  112. #        local i,v;
  113. #           v := pv(p);
  114. #           for i in v do
  115. #              if op(1,i) = 'v1(t)' then RETURN(op(2,i)) fi
  116. #           od;
  117. #           ERROR(`Could not find v1`);
  118. #      end;
  119. #-----------------------------
  120. #
  121. # More to be recommended is the "trick" of substituting the solution
  122. # (a set of equations) into the desired variable name.
  123. #
  124. > V1 := T -> subs(pv(T), v1(t));
  125.                          V1 := T -> subs(pv(T), v1(t))
  126.  
  127. > V1(1);
  128.                                   1.599162730
  129.  
  130. > V1(1.5);
  131.                                   1.579211333
  132.  
  133. > quit
  134.