gnuplot-beta

integration approximation demo in bivariat.dem

[Daniel J S. <dan...@ie...>]

Paging through all.dem, something caught my eye in the bivariat.dem demo.  It's not of critical importance; after all, gnuplot is about plotting not applications.

The first example is of approximating integration of a function.  Unfortunately it appears to be a bad approximation.  Attached is a PNG of the plot.  Note how the integral approximation strangely drifts downward as x tends to +5.  This is an integral of a strictly positive function!

Anyway, given the formula for approximation (i.e., via rectangles), one has to conclude it is a sampling sort of thing.

So

1) I changed the number of samples so that x=0 ends up being one of the samples.  That corrected the slow drift away from the asymptote.

2) ... but introduced a problem where gnuplot's choice of samples and the x<=0 test caused a strange artifact.  To fix that, I used x<=epsilon.

3) ... and there still remained an offset of delta/2 due to the integration process, so I added in the appropriate offset.

See the attached PNG and diff files for an idea of the needed changes.

Anyway, do we want to fix these sorts of things or just live with it?

Dan

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

Re: integration approximation demo in bivariat.dem

[Daniel J S. <dan...@ie...>]

Ah!  Very nice.  Thank you.  This should be moved into CVS as soon as possible.  Developers are probably being overwhelmed with patches, so perhaps place in SourceForge if it isn't attended to in the near future.

Dan


James R. Van Zandt wrote:
> 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.

Not following what you mean.  Do you mean use a lower order approximation such as trapezoidal to compare against Simpson's rule?


>                - 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.

Very nice.  Could add a note about the accuracy of the third order Simpson's for approximating a second order function at the sample points, i.e., exact.  Is that the kind of thing you meant?


>  - 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.

I didn't go much beyond the first, but studying it, that jumps out now.


>  - For the last plot, use "with points" to make the patterns more
> apparent, and because the function is only defined over the integers.

Much better.

Dan

Re: integration approximation demo in bivariat.dem

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

Daniel J Sebald <dan...@ie...> wrote:
>   James R. Van Zandt wrote:
>   > 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.
>
>   Not following what you mean.  Do you mean use a lower order
>   approximation such as trapezoidal to compare against Simpson's
>   rule?

I'm just thinking it would be easier for the user to study the
original program first, then the more sophisticated one.
...
>   Very nice.  Could add a note about the accuracy of the third order
>   Simpson's for approximating a second order function at the sample
>   points, i.e., exact.  Is that the kind of thing you meant?

Yes.

             - Jim Van Zandt