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