1
2critical_points(expr,var):= block([numer:true,f:diff(expr,var)],
3    realroots(num(f)*denom(f),10^-7));
4
5function_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);
11upper_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
19graph_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
51upper_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
56make_partition(a,b,n):=makelist(a+(b-a)*i/n,i,0,n);
57
58