On 27 Aug 2008, Mahery Raharinjatovo wrote:
>I would want to plot the curve define by points of coordinates (t,v)
>
>chute():=block([g,k,r,rho,rhoh,vm,V,v],
>g:9.8,
>k:9.33/100,
>rho:2400,
>rhoh:950,
>r:1.5/1000,
>V:(4/3)*3.14*r*r*r,
>vm:(g*rho*V/k)*(1-rhoh/rho),
>dt:1.5/10000,
>a:-k/(rho*V),
>b:g*(1-(rhoh/rho)),
>v:0,t:0,
>while v<0.999*vm do (print(v,"?",t),dv:(a*v+b)*dt,v:v+dv,t:t+dt)
>)$
>(%i20) chute();
>0 ? 0
>8.8812499999999992E-4 ? 1.4999999999999999E-4
>0.0014097334461253 ? 2.9999999999999997E-4
>etc, etc.
Here is a "plain jane" solution similar to traditional
procedural programming. (There are fancier ways
to do almost everything using Maxima!).
1. create a batch file which I called pnts1.mac.
Here is the file:
===============
/* file pnts1.mac */
( fpprintprec : 8, display2d : false)$
(g : 9.8, k : 9.33e-2, rho : 2400.0, rhoh : 950.0, r : 1.5e-3,
vcap : float(4*%pi*r^3/3), vm : (g*rho*vcap/k)*(1-rhoh/rho),
b : g*(1-(rhoh/rho)), a : -k/(rho*vcap), vmax : 0.999*vm )$
print( g, k, rho, rhoh)$
print( r, vcap, vm, b)$
print( a, vmax)$
/* now let's make the needed list */
(d t : 1.5e-4, t : 0, v : 0, tlist : [0], vlist : [0],
print( t, v),
while ( v < vmax ) do (
t : t + dt, dv : dt*(a*v + b), v : v + dv,
print( t, v),
tlist : cons( t, tlist),
vlist : cons( v, vlist) ),
tvlist : makelist( [ tlist[i], vlist[i] ], i, 1, length(tlist) ),
tvlist : reverse( tvlist) )$
tvlist;
( load(draw), load(qdraw),
qdraw( xr(-2e-4, 2.5e-3 ), yr( -2e-4, 2.5e-3 ),
pts( tvlist, ps(2), pc(blue) ) ) )$
====================
Here is output if I use batch(pnts1)$
========================
(%i1) batch(pnts1)$
batching #pc:/work3/pnts1.mac
(%i2) (fpprintprec : 8, display2d : false)
(%i3)
(g:9.8,k:0.0933,rho:2400.0,rhoh:950.0,r:0.0015,vcap:float(4*%pi*r^3/3),
vm:g*rho*vcap*(1-rhoh/rho)/k,b:g*(1-rhoh/rho),a:(-k)/(rho*vcap),
vmax:0.999*vm)
(%i4) print(g,k,rho,rhoh)
9.8 0.0933 2400.0 950.0
(%i5) print(r,vcap,vm,b)
0.0015 1.41371669E-8 0.00215315 5.9208333
(%i6) print(a,vmax)
-2749.8437 0.002151
(%i7) (dt:1.5E-4,t:0,v:0,tlist:[0],vlist:[0],print(t,v),
while v < vmax do
(t:dt+t,dv:dt*(b+a*v),v:dv+v,print(t,v),tlist:cons(t,tlist),
vlist:cons(v,vlist)),
tvlist:makelist([tlist[i],vlist[i]],i,1,length(tlist)),
tvlist:reverse(tvlist))
0 0
1.5E-4 8.88125E-4
3.0E-4 0.00140992
4.5E-4 0.00171649
6.0E-4 0.0018966
7.5E-4 0.00200242
9.0E-4 0.0020646
0.00105 0.00210112
0.0012 0.00212258
0.00135 0.00213519
0.0015 0.0021426
0.00165 0.00214695
0.0018 0.00214951
0.00195 0.00215101
(%i8) tvlist
(%o8) [[0,0],[1.5E-4,8.88125E-4],[3.0E-4,0.00140992],[4.5E-4,0.00171649],
[6.0E-4,0.0018966],[7.5E-4,0.00200242],[9.0E-4,0.0020646],
[0.00105,0.00210112],[0.0012,0.00212258],[0.00135,0.00213519],
[0.0015,0.0021426],[0.00165,0.00214695],[0.0018,0.00214951],
[0.00195,0.00215101]]
(%i9) (load(draw),load(qdraw),
qdraw(xr(-2.0E-4,0.0025),yr(-2.0E-4,0.0025),
pts(tvlist,ps(2),pc(blue))))
" qdraw(...), qdensity(...), syntax: type qdraw()\; "
/* plot appears in separate window */
(%i11)
==================
the requested plot comes up immediately
as the last action of this batch file.
-------------------------------------------
For the plot, I used my qdraw package (available on
my webpage with chapter 5 of Maxima
by Example) which is a simple interface to
draw2d. Naturally, you can use your own
favorite plotting method.
----------------
Note also that you need "float(...)" of an
expression involving %pi to get the
%pi converting to floating point.
===============
Before looking at the above code,
try out (interactively) the following prototype
code for creating a list of
(x, w) pairs:
================
(%i44) (w : 0.0, x : 0.0, xlist : [0], wlist : [0],
while float(w) < 20.0 do(
x:x + 1.0,
w : w+x^2, print( x, w),
xlist : cons( x, xlist),
wlist : cons( w, wlist) ),
xwlist : makelist( [ xlist[i], wlist[i] ], i, 1, length(xlist) ),
xwlist : reverse( xwlist) );
1.0 1.0
2.0 5.0
3.0 14.0
4.0 30.0
(%o44) [ [0, 0], [1.0, 1.0], [2.0, 5.0], [3.0, 14.0], [4.0, 30.0]]
==================
Ted Woollett
http://www.csulb.edu/~woollett