Menu
▾
▴
gnuplot-beta
|
From: 3snoW <vas...@gm...> - 2011-06-16 00:36:00
|
Hello, I'm not sure if this qualifies as "Dev", If it doesn't I'm sorry. I've searched the web for ways to integrate in gnuplot and only found people saying it is not possible. Unsatisfied, I decided to see if i could make one myself. I could :D This was the result: 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); Basically, intfAux 1 and 2 are just functions for this thing to work, and intf(a,b) integrates f(x) from a to b. You can just type this in every time you want to integrate f(x), or you can download the gnuplot.ini that i made with this function and put it into your binary folder: http://old.nabble.com/file/p31856137/gnuplot.ini gnuplot.ini Here is also a demonstration of using this: http://old.nabble.com/file/p31856137/Integrate%2B.png There are two problems with this function (that I know of): -It will most likely not work with self recursive functions because this function itself is self recursive, and gnuplot has a "stack overflow" limit. This makes it impossible to double integrate using this function. I made it so it was not at the limit, it is almost there, so, for example having a function g(x)=h(intf(0,x)) should not be a problem. By the way, if anyone knows of a way to set the stack overflow limit higher, I'd appreciate posting it here. -The other problem is that it will integrate specifically f(x), so if you want to integrate another function, g(x) for example, you would have to create an intg(a,b). To make this task easier, I made it so that if you replace in those 3 lines of code all the f's with the name of your function, g for example (again), you would have your intg made. This way you can just open notepad, use the "replace" command and copy the result to gnuplot. 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! 3snoW -- View this message in context: http://old.nabble.com/Integration-in-gnuplot-IS-possible%21-tp31856137p31856137.html Sent from the Gnuplot - Dev mailing list archive at Nabble.com. |
|
From: Daniel J S. <dan...@ie...> - 2011-06-16 08:37:43
|
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 |
|
From: 3snoW <vas...@gm...> - 2011-06-16 10:54:25
|
Daniel J Sebald wrote: > > 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. > > [snip] > > 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 > Hi Dan, The reason I was using f(x)-f(a) in the intfAux functions and then add f(a)*(x-a) is because I originally was using a fixed small step h instead of dividing the interval in 90 small parts. That way, if f(a) was not 0, the integral would have a discontinuity every time x was a multiple of h, so, as a workaround to this problem, I set it to integrate f(x)-f(a), and compensated adding f(a)*(x-a). But since I switched to dividing into 90 equal parts, this is no longer a problem. I didn't make the integration method more sophisticated because of the 'stack overflow' issue. Adaptive integral techniques would most likely consume 'stacks' at a much higher rate, and in the end you would have a worse integration because of the lack of resolution. I'm thinking of maybe upgrading this to the trapezoidal integration or even a superior order technique, but for now I think this works fine for any 'well behaved' function. Lastly, the generalized function that includes N as a parameter is a good idea, but I'll call it Nintf instead of intf, because that way I can have both functions. Thank you for your contribution, it was really constructive! I'll update the config.ini file with your changes. 3snoW -- View this message in context: http://old.nabble.com/Integration-in-gnuplot-IS-possible%21-tp31856137p31859432.html Sent from the Gnuplot - Dev mailing list archive at Nabble.com. |
|
From: Hans-Bernhard B. <HBB...@t-...> - 2011-06-16 18:56:31
|
On 16.06.2011 02:35, 3snoW wrote: > I'm not sure if this qualifies as "Dev", If it doesn't I'm sorry. > I've searched the web for ways to integrate in gnuplot and only found people > saying it is not possible. Strange --- particularly since one of our published/enclosed demos, bivariat.dem, has been doing just that since 1998! |
|
From: Daniel J S. <dan...@ie...> - 2011-06-16 20:06:38
|
On 06/16/2011 01:57 PM, Hans-Bernhard Bröker wrote: > On 16.06.2011 02:35, 3snoW wrote: > >> I'm not sure if this qualifies as "Dev", If it doesn't I'm sorry. >> I've searched the web for ways to integrate in gnuplot and only found people >> saying it is not possible. > > Strange --- particularly since one of our published/enclosed demos, > bivariat.dem, has been doing just that since 1998! I gave those formulas a try. They behave better than 3snoW's formulas in the example I used, because of the parameter delta (i.e., h) rather than the fixed number of intervals of 3snoW and because of its use of Simpson's rule. The bivariat.dem formula still has the problem of recursion/stack depth, in this case: recursion depth limit exceeded An internal integration function is probably the only solution to that. The tricky part would be interpreting the function name, i.e., f(x) = cos(x) h = 0.2 M = 1 plot intgrt(f(x),h,M) where f(x) is the function, h is the (initial) stepsize, M is the integration method. The first argument, in this case "f(x)", would need special treatment in the parser different from other functions. Dan |
|
From: Ethan A M. <sf...@us...> - 2011-06-16 21:18:22
|
On Thursday, June 16, 2011 01:06:23 pm Daniel J Sebald wrote: > The bivariat.dem formula still has the problem of > recursion/stack depth, in this case: > > recursion depth limit exceeded The depth limit was set rather arbitrary at 250 on the logic that "nobody would need more than that, right?" It would be easy enough to make it user-defined. If you increase it to something much larger and the program runs out of stack space, then at least you know it's your own fault. Ethan |
|
From: Daniel J S. <dan...@ie...> - 2011-06-19 07:58:28
|
On 06/16/2011 03:59 PM, Ethan A Merritt wrote: > On Thursday, June 16, 2011 01:06:23 pm Daniel J Sebald wrote: >> The bivariat.dem formula still has the problem of >> recursion/stack depth, in this case: >> >> recursion depth limit exceeded > > The depth limit was set rather arbitrary at 250 on the logic that > "nobody would need more than that, right?" > > It would be easy enough to make it user-defined. > If you increase it to something much larger and the program > runs out of stack space, then at least you know it's your own > fault. I suppose that is a solution. Or could remove the limit and see how gnuplot behaves when memory is used up. That solution doesn't address the inefficient method of integration, but seeing as the feature probably isn't used often the solution might work well enough. Dan |
|
From: Pieter-Tjerk de B. <ptd...@cs...> - 2011-06-16 23:03:03
|
Hello, On Wed, Jun 15, 2011 at 05:35:54PM -0700, 3snoW wrote: > I've searched the web for ways to integrate in gnuplot and only found people > saying it is not possible. Unsatisfied, I decided to see if i could make one > myself. I could :D Nice! > Let me know if you have any ideas to make it better I'd suggest using the recursion not to cover the interval _linearly_, but using the recursion to _bisect_ the interval, like this: intf(a,x,n) = (n==0) ? f((a+x)/2)*(x-a) : ( intf(a,(a+x)/2,n-1) + intf((a+x)/2,x,n-1) ) The last parameter n is the number of recursions, and which equals the 2log of the number of points in which to cut the interval. This approach should practically avoid stack overflows, since much less recursion depth is needed for the same accuracy. Regards, Pieter-Tjerk |
|
From: Alexandre F. <o.a...@gm...> - 2011-06-17 11:30:50
|
Does someone read what i have wrote? |
|
From: Alexandre F. <o.a...@gm...> - 2011-06-17 13:01:08
|
The Pieter-Tjerk is fairly better :D redefining it this make the same thing but you ask for a maximum integration distance instead a maximum intf(x0,a,h) = (h<=delta) ? f(x0+a/2)*a : ( intf(x0,a/2,delta) + intf(x0+a/2,a/2,delta) ) when you use gnuplot i think you want to plot something, this version if you say plot intf(x0, x - x0, delta) ensures that any integration interval will be larger than delta, and that the recurrence steep will be the minimum (with this function structure) if you plot with 1024 samples the integration function intf(x0, a, delta) is called 1397077 times, and the function f(x) will be evaluated 350548 times and the max recurrence depth will be 10 to integrate with the same acuracy the Daniel's scheme will call the intf* functions and f(x) are called about 523 thounsand times. and the max recurrence deepth is 1024 my scheme calls f(x) 1024 times and intf(x) 1024 times, but use some global variables. excuse-me i accidentally replied only to 3snoW here goes the my scheme. It's not encouraged the use of global variables, but in such case i think here we a better solution, not THE better solution # Characteristic for the integration n_samps = 100; a = 0.0; b = 5.0; # initialization aux = 0.0; h = (b-a)/n_samps; set samples n_samps set xrange [a to b] # the function f(x) = sin(5*x)/(0.1+x) #the integrator scheme intf(x) = (aux = (aux + f(x)*h)) # plot it plot f(x), intf(x) # another function f(x) = sin(x)*x replot for the 2D is also possible. the only thing you need is to reinit the aux variable when x=a this scheme is as poor as the former, therefore it's O(nsamples) and why to use recurrence if you can use memorization? intf can also have other integration schemes, with the advantage that you have to implement only once for example simpson's rule or any other method intf(x) = (aux = (aux + (f(x-h)+4*f(x-0.5*h)+f(x))*h/ 6)) this example make the h function to be called three times more. the ideal to generate plots is to derive integration schemes open in f(x-h) and closed in f(x), but i will not make this now :D |
|
From: Alexandre F. <o.a...@gm...> - 2011-06-17 13:01:56
|
The Pieter-Tjerk is fairly better :D redefining it this make the same thing but you ask for a maximum integration distance instead a maximum intf(x0,a,h) = (h<=delta) ? f(x0+a/2)*a : ( intf(x0,a/2,delta) + intf(x0+a/2,a/2,delta) ) when you use gnuplot i think you want to plot something, this version if you say plot intf(x0, x - x0, delta) ensures that any integration interval will be larger than delta, and that the recurrence steep will be the minimum (with this function structure) if you plot with 1024 samples the integration function intf(x0, a, delta) is called 1397077 times, and the function f(x) will be evaluated 350548 times and the max recurrence depth will be 10 to integrate with the same acuracy the Daniel's scheme will call the intf* functions and f(x) are called about 523 thounsand times. and the max recurrence deepth is 1024 my scheme calls f(x) 1024 times and intf(x) 1024 times, but use some global variables. excuse-me i accidentally replied only to 3snoW here goes the my scheme. It's not encouraged the use of global variables, but in such case i think here we a better solution, not THE better solution # Characteristic for the integration n_samps = 100; a = 0.0; b = 5.0; # initialization aux = 0.0; h = (b-a)/n_samps; set samples n_samps set xrange [a to b] # the function f(x) = sin(5*x)/(0.1+x) #the integrator scheme intf(x) = (aux = (aux + f(x)*h)) # plot it plot f(x), intf(x) # another function f(x) = sin(x)*x replot for the 2D is also possible. the only thing you need is to reinit the aux variable when x=a this scheme is as poor as the former, therefore it's O(nsamples) and why to use recurrence if you can use memorization? intf can also have other integration schemes, with the advantage that you have to implement only once for example simpson's rule or any other method intf(x) = (aux = (aux + (f(x-h)+4*f(x-0.5*h)+f(x))*h/6)) this example make the h function to be called three times more. the ideal to generate plots is to derive integration schemes open in f(x-h) and closed in f(x), but i will not make this now :D |
|
From: Alexandre F. <o.a...@gm...> - 2011-06-17 13:03:43
|
Pieter-Tjerk method is fairly better :D redefining it this make the same thing but you ask for a maximum integration distance instead a maximum intf(x0,a,h) = (h<=delta) ? f(x0+a/2)*a : ( intf(x0,a/2,delta) + intf(x0+a/2,a/2,delta) ) when you use gnuplot i think you want to plot something, this version if you say plot intf(x0, x - x0, delta) ensures that any integration interval will be larger than delta, and that the recurrence steep will be the minimum (with this function structure) if you plot with 1024 samples the integration function intf(x0, a, delta) is called 1397077 times, and the function f(x) will be evaluated 350548 times and the max recurrence depth will be 10 to integrate with the same acuracy the Daniel's scheme will call the intf* functions and f(x) are called about 523 thounsand times. and the max recurrence deepth is 1024 my scheme calls f(x) 1024 times and intf(x) 1024 times, but use some global variables. excuse-me i accidentally replied only to 3snoW here goes the my scheme. It's not encouraged the use of global variables, but in such case i think here we a better solution, not THE better solution # Characteristic for the integration n_samps = 100; a = 0.0; b = 5.0; # initialization aux = 0.0; h = (b-a)/n_samps; set samples n_samps set xrange [a to b] # the function f(x) = sin(5*x)/(0.1+x) #the integrator scheme intf(x) = (aux = (aux + f(x)*h)) # plot it plot f(x), intf(x) # another function f(x) = sin(x)*x replot for the 2D is also possible. the only thing you need is to reinit the aux variable when x=a this scheme is as poor as the former, therefore it's O(nsamples) and why to use recurrence if you can use memorization? intf can also have other integration schemes, with the advantage that you have to implement only once for example simpson's rule or any other method intf(x) = (aux = (aux + (f(x-h)+4*f(x-0.5*h)+f(x))*h/6)) this example make the h function to be called three times more. the ideal to generate plots is to derive integration schemes open in f(x-h) and closed in f(x), but i will not make this now :D |
