Re: Integration in gnuplot IS possible!

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

On 06/15/2011 07:35 PM, 3snoW wrote:
>
> Hello,
> I'm not sure if this qualifies as "Dev", If it doesn't I'm sorry.
[snip]
> Let me know if you have any ideas to make it better or if you find any bug
> using it! I hope this was helpful!

Very creative.  Of course, this is straightforward rectangular 
approximation to integration.  It will have problems with functions that 
change quickly near discontinuities and so on so won't have the 
precision of, say, adaptive integral techniques.  Nonetheless, it works 
(sort of, see note at end), assuming the initial condition is that the 
integral passes through zero at a.

Let's see if we can simplify your formulas slightly for efficiency. 
Starting with

intfAux1(a,x,h)=(x>a?(f(x)-f(a))*h+intfAux1(a,x-h,h):0);
intfAux2(a,x,h)=(x<a?(f(x)-f(a))*h+intfAux2(a,x+h,h):0);
intf(a,x)=x>a?intfAux1(a,x-abs(x-a)/90./2,abs(x-a)/90.)+f(a)*(x-a):-intfAux2(a,x+abs(x-a)/90./2,abs(x-a)/90.)+f(a)*(x-a);

Note that the abs() function call is extraneous because we know that

If x>a, then abs(x-a) = x-a
If x<=a, then abs(x-a) = a-x

Therefore, substituting the above

intf(a,x)=x>a?intfAux1(a,x-(x-a)/90./2,(x-a)/90.)+f(a)*(x-a):-intfAux2(a,x+(a-x)/90./2,(a-x)/90.)+f(a)*(x-a);

It's possible to factor out the h from the auxiliary formulas and simply 
multiply it with the auxiliary formula, knowing that h=(x-a)/90 in the 
one case and h=(a-x)/90 in the other case:

intfAux1(a,x,h)=(x>a?(f(x)-f(a))+intfAux1(a,x-h,h):0);
intfAux2(a,x,h)=(x<a?(f(x)-f(a))+intfAux2(a,x+h,h):0);
intf(a,x)=x>a?intfAux1(a,x-(x-a)/90./2,(x-a)/90.)*(x-a)/90.+f(a)*(x-a):-intfAux2(a,x+(a-x)/90./2,(a-x)/90.)*(a-x)/90.+f(a)*(x-a);

Swap the order of x and a in the second formula:

intf(a,x)=x>a?intfAux1(a,x-(x-a)/90./2,(x-a)/90.)*(x-a)/90.+f(a)*(x-a):intfAux2(a,x+(a-x)/90./2,(a-x)/90.)*(x-a)/90.+f(a)*(x-a);

Note that f(a) is subtracted 90 times, so that can be removed from the 
auxiliary formulas and lumped in the main formula:

intfAux1(a,x,h)=(x>a?f(x)+intfAux1(a,x-h,h):0);
intfAux2(a,x,h)=(x<a?f(x)+intfAux2(a,x+h,h):0);
intf(a,x)=x>a?intfAux1(a,x-(x-a)/90./2,(x-a)/90.)*(x-a)/90.-90.*f(a)*(x-a)/90.+f(a)*(x-a):intfAux2(a,x+(a-x)/90./2,(a-x)/90.)*(x-a)/90.-90.*f(a)*(x-a)/90.+f(a)*(x-a);

at which point a couple terms cancel.  So we have

intfAux1(a,x,h)=(x>a?f(x)+intfAux1(a,x-h,h):0);
intfAux2(a,x,h)=(x<a?f(x)+intfAux2(a,x+h,h):0);
intf(a,x)=x>a?intfAux1(a,x-(x-a)/90./2,(x-a)/90.)*(x-a)/90.:intfAux2(a,x+(a-x)/90./2,(a-x)/90.)*(x-a)/90.;

Seems to work the same.  It's not drastically shorter, but the 
arithmetic operations are reduced.

The thing one has to be aware of is that as x gets further and further 
away from a, the resolution effectively decreases because element 
h=(x-a)/90 becomes larger.  This is a problem for functions that appear 
to be high frequency.  For example, try your cosine example, but with

set sample 1000
set xrange [-200:200]

Notice how accuracy is lost moving outward.  So you might want to 
experiment with making h another parameter in the function.  Generalizing,

intf(a,x,N)=x>a?intfAux1(a,x-(x-a)/N/2,(x-a)/N)*(x-a)/N:intfAux2(a,x+(a-x)/N/2,(a-x)/N)*(x-a)/N;

I can't push N high enough to get good resolution before a stack 
overflow occurs.

An internal integration feature is probably the only good way to address 
this sort of thing.

Dan