share/numeric/riemsum.mac

   1
   2 critical_points(expr,var):= block([numer:true,f:diff(expr,var)],
   3     realroots(num(f)*denom(f),10^-7));
   4
   5 function_max(expr,var,a,b,critical_points):=
   6    block([y:max(at(expr,var=a),at(expr,var=b)) ],
   7           for v in critical_points
   8              do (v:at(x,v),
   9              if a<v and v<b and at(expr,var=v) > y then y:at(expr,var=v)),
  10           y);
  11 upper_sum(expr,var,partition):=
  12   block([crit:critical_points(expr,var)],
  13     partition:sort(partition),
  14      sum(function_max(expr,var,partition[i],partition[i+1],crit)*
  15              (partition[i+1]-partition[i]),i,1,length(partition)-1));
  16
  17 /* this will be part of maxima in a little while .. */
  18
  19 graph_riemann_sum(up,expr,var,partition):=
  20   block([lis,numer:true,crit:critical_points(expr,var),n:length(partition),xx,
  21          interval,display2d:false,sgn,min:100000,max:-10000000],
  22     sgn: if up= upper then 1 else -1,
  23     partition:sort(partition),
  24     interval:partition[n]-partition[1],
  25     lis:[],
  26     for i:0 thru 50 do (xx:partition[1]+i*interval/50,
  27                         lis:cons(val:at(expr,var=xx),lis),
  28                         lis:cons(xx,lis),
  29                         if val < min then min:val,
  30                         if val > max then max:val),
  31     sprint("{plot2d ", "{xrange " ,partition[1]-.5,partition[n]+.5,"} {yrange" ,
  32            min-.5,max+.5, "} {xversusy "),
  33     tcl_output(lis,1),
  34     tcl_output(lis,2),
  35     sprint( "}"),
  36     sprint(" {xversusy {"),
  37     for i:1 thru n-1
  38      do(
  39        sprint(partition[i],partition[i],partition[i+1],partition[i+1],
  40               partition[i])),
  41     sprint(" } { "),
  42     for i:1 thru n-1
  43       do(xx:sgn*function_max(sgn*expr,var,partition[i],partition[i+1],crit),
  44          sprint(0.0,xx,xx,0.0,0.0)),
  45     sprint("}}}"),
  46     ""
  47       );
  48
  49
  50
  51 upper_and_lower_sums(expr,var,partition):=
  52    block([up:upper_sum(expr,var,partition),
  53           low:-upper_sum(-expr,var,partition)],
  54     [up,low,up-low]);
  55
  56 make_partition(a,b,n):=makelist(a+(b-a)*i/n,i,0,n);
  57