Re: integration approximation demo in bivariat.dem

[James R. V. Z. <jr...@co...>]

Daniel J Sebald <dan...@ie...> wrote:
>   Paging through all.dem, something caught my eye in the bivariat.dem demo.  
>   ...
>   The first example is of approximating integration of a function.
>   Unfortunately it appears to be a bad approximation.

I looked the integration demo some time ago.  Attached are some
suggested changes.  It might be worth while including two versions -
one to show the basic approach, and the other to show possible
refinements that improve the results at the cost of some obscurity.

               - Jim Van Zandt

These are the changes:
 - Adjusting the step size so the range of integration is an integral
number of steps.
 - Using Simpson's rule rather than just taking the function value at
one end of each step.
 - Define the first integrand immediately before the plot command.
 - Add points from the real erf(x) to the first plot, for comparison.
 - For the second plot, define the integrand f(x)=cos(x).  That way,
the integral is sin(x) instead of being offset by a constant of
integration.
 - For the last plot, use "with points" to make the patterns more
apparent, and because the function is only defined over the integers.
 - Update the "set key" command to the current syntax


--- demo/bivariat.dem.orig	2006-03-21 21:33:02.000000000 -0500
+++ demo/bivariat.dem	2006-06-21 19:57:40.000000000 -0400
@@ -9,35 +9,42 @@
 # integral2_f(x,y) approximates the integral from x to y.
 # define f(x) to be any single variable function
 #
-# the integral is calculated as the sum of f(x_n)*delta 
-#   do this x/delta times (from x down to 0)
+# the integral is calculated using Simpson's rule as 
+#          ( f(x-delta) + 4*f(x-delta/2) + f(x) )*delta/6
+# repeated x/delta times (from x down to 0)
 #
-f(x) = exp(-x**2)
 delta = 0.2
 #  delta can be set to 0.025 for non-MSDOS machines
 #
 # integral_f(x) takes one variable, the upper limit.  0 is the lower limit.
 # calculate the integral of function f(t) from 0 to x
-integral_f(x) = (x>0)?integral1a(x):-integral1b(x)
-integral1a(x) = (x<=0)?0:(integral1a(x-delta)+delta*f(x))
-integral1b(x) = (x>=0)?0:(integral1b(x+delta)+delta*f(x))
+# choose a step size no larger than delta such that an integral number of
+# steps will cover the range of integration.
+integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta))
+int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.)
+int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.)
 #
 # integral2_f(x,y) takes two variables; x is the lower limit, and y the upper.
-# claculate the integral of function f(t) from x to y
-integral2_f(x,y) = (x<y)?integral2(x,y):-integral2(y,x)
-integral2(x,y) = (x>y)?0:(integral2(x+delta,y)+delta*f(x))
+# calculate the integral of function f(t) from x to y
+integral2_f(x,y) = (x<y)?int2(x,y,(y-x)/ceil((y-x)/delta)): \
+                        -int2(y,x,(x-y)/ceil((x-y)/delta))
+int2(x,y,d) = (x>y-d*.5) ? 0 : (int2(x+d,y,d) + (f(x)+4*f(x+d*.5)+f(x+d))*d/6.)
 
 set autoscale
 set title "approximate the integral of functions"
 set samples 50
 
-plot [-5:5] f(x) title "f(x)=exp(-x**2)", 2/sqrt(pi)*integral_f(x) title "erf(x)=2/sqrt(pi)*integral_f(x)"
+f(x) = exp(-x**2)
+
+plot [-5:5] f(x) title "f(x)=exp(-x**2)", \
+  2/sqrt(pi)*integral_f(x) title "erf(x)=2/sqrt(pi)*integral_f(x)", \
+  erf(x) with points
 
 pause -1 "Hit return to continue"
 
-f(x)=sin(x)
+f(x)=cos(x)
 
-plot [-5:5] f(x) title "f(x)=sin(x)", integral_f(x)
+plot [-5:5] f(x) title "f(x)=cos(x)", integral_f(x)
 
 pause -1 "Hit return to continue"
 
@@ -78,7 +85,7 @@
 set yrange [-10:10]
 set isosamples 10
 set samples 100
-set key at 4,-3
+set key 4,-3
 set title "Min(x,y) and Max(x,y)"
 
 #
@@ -104,7 +111,7 @@
 
 set title "Greatest Common Divisor (for integers only)"
 
-plot gcd(x, 60)
+plot gcd(x, 60) with points
 pause -1 "Hit return to continue"
 
 reset