High Accuracy Technique For Solving Numerical Integrals

ERROR AREA REDUCTION INTEGRAL 

In this material, I want to show a way to compute the definite integral of a function y = f(x), between x = x1 and x = x2 with accuracy greater than both Simpson's rule and the polynomial rule.

Consider a polynomial approximation for the curve y=f(x), say y = h(x). where h(x) is a polynomial function that approximates f(x).

Since f(x) is the function which we wish to integrate and so is already known, we need to compute h(x) which is our ?

  Approximating function.

Sure you knew that already!

I trust that there are countless methods for generating the polynomial function h(x), but here we shall use one of the simplest ones available.

First we choose the limits for the approxiation say x = x1 and x = x2  between which we wish to compute the integral

Next, we choose the highest integer power that we wish to have in the polynomial, say n, which will be the degree of our polynomial, h(x).
So the polynomial will be h(x) = a₀ + a₁x + a₂x² + a₃x³ ....+ a_n-1x^n-1 + a_nxⁿ

Now, we can easily calculate the values of a₀ , a₁ a₂ a_{3 .....}a_n 
since we know the value of the function we wish to approximate.

All we do is to take a number of collision steps, between x = x₁ and x = x₂   , say m steps.

These steps define the exact x-collision values where y=f(x) touches its approximating polynomial,
y=h(x)

m is directly related to the degree of our polynomial approximation by m = n+1

Now we calculate the x-spacing between each of the steps, say dx.

The x points we are taking are:

x1, x1+dx, x1+2.dx, x1+3.dx, ....+ x1+n.dx.
So since x₂   = x₁ + n.dx.

Then, dx = (x2 - x1)/n.

So if we wish to approximate y=sin(x), for example: say with a polynomial of degree n = 5,

between x = 1 and x = 3.

Then we would have:

dx = (3 - 1)/5; say dx = 0.4;

So x1 = 1, x2 = 1.4, x3 = 1.8, x4 = 2.2, x5 = 2.6, x6 = 3.

Of course this is consistent with the number of collision points, e.g m being n+1 = 6.

So the approximating polynomial gives:

 h(1.0) = sin(1)   = a0 + 1.0.a1 + 1.02.a2 + 1.03.a3 + 1.04.a4 + 1.05.a5

 h(1.4) = sin(1.4) = a0 + 1.4.a1 + 1.42.a2 + 1.43.a3 + 1.44.a4 + 1.45.a5

 h(1.8) = sin(1.8) = a0 + 1.8.a1 + 1.82.a2 + 1.83.a3 + 1.84.a4 + 1.85.a5

 h(2.2) = sin(2.2) = a0 + 2.2.a1 + 2.22.a2 + 2.23.a3 + 2.24.a4 + 2.25.a5

 h(2.6) = sin(2.6) = a0 + 2.6.a1 + 2.62.a2 + 2.63.a3 + 2.64.a4 + 2.65.a5

 h(3.0) = sin(3.0) = a0 + 3.0.a1 + 3.02.a2 + 3.03.a3 + 3.04.a4 + 3.05.a5

You see where this is going?

                                     Sure I do!

This is a perfect system of linear equations which can be solved very quickly by a Matrix program!
And guess what? I  have one for you if you are really interested! and what's more?

 its absolutely free! 

I have programmed them in the Java and C languages and once requests are made in the comments section, I will make the code available.

The system's solution is then:

h(x) = 0.05621527890686895 + 0.8070522394224585.x + 0.2659311952212848.x2 - 0.35406795117551315.x3 + 0.0696942723181288.x4 - 0.003354049885331352.x5

Try and find say h(1.7), bearing in mind that 1.7 lies between  x = 1 and x = 3, the validity limits for the approximation.

You get h(1.7) = 0.9916802157371398

while its real value, which is f(1.7) = sin(1.7) = 0.9916648104524686

Need I tell that if we pick a larger degree for n, we will get a longer and more accurate approximating function, h(x) for

For instance for n = 10, we get

h(1.7) = 0.9916648104538464

which is much closer to:

f(1.7) = sin(1.7) = 0.9916648104524686


Now the gist of all this is that since h(x) approximates f(x),. we can apply simple polynomial integration to integrate h(x) since

h(x) ≅ f(x)




So also:

∫h(x)dx ≅ ∫f(x)dx


Now for a bit of theory, just a light dose, nothing too complicated!

Remember that diagram? the one about a naughty curve up there? well here it is!




The area under the whole curve, i.e from A_{1 to}  A₂ along y=f(x) is the definite integral of   y=f(x) between  x1 and  x2.

The area between A₁ and A₂  along y=h(x) is the definite integral of the approximating function which is hence an approximate value for the definite integral of y=f(x).

It is evident then that the area A_e between y=f(x) and y=h(x) is the error incurred in computing
the definite integral of y=f(x) with the polynomial integral of y=h(x).




The point of discovery!

So can we compute or estimate the error area, A_e ?

If we can compute it, then by adding it to A_h , we can have a better approximation for the definite integral of  y=f(x).

If the absolute definite integral is A_f , then 

A_f  = A_h + A_e

Which can also be represented as:

Before we proceed, a very important principle to note is that any approximating function can be applied to compute h(x). Also any ingenious thinker can come up with a way to compute A_e

I propose a method for calculating A_e here in this material.

Our technique will involve Simpson's rule.

First, we use Simpson's rule to estimate A_f  .

Simpson's approximation for ∫f(x)dx between x1 and  x2 is:

Next, we use Simpson's rule to estimate A_h . Simpson's approximation for ∫h(x)dx is:

Can you see where this is going? Maybe not exactly this time?

Well, from

A_f  = A_h + A_e .............(1)

This gives us that:

A_e  = A_f  - A_h

Now, since S(f(x)) is the Simpson's area approximation for A_{f and} S(h(x)) is the Simpson's area approximation for A_h , let S_{f =} S(f(x)) and S_h = S(h(x)).

Then, we see from the graph that:

S_f  = S_h + S_e.............(2)

S_e  = S_f - S_h.............(3)

where S_e is Simpson's approximation for the error area, A_e.

Now, we assume that since A_e ≅ S_e, so we can substitute equation (3) to replace S_e in equation (1).

Then we have:

A_f = A_h + S_f - S_{h.........(4)}

And we are nigh done!

For A_f =

and A_h =

and S_f - S_h is equal to the difference between:

AND 

This gives S_f - S_h =

Surely you can see the silver lining behind that cloud now?


Yes, but for one last step!

Now for a larger area comprising many of the above areas, we take the sum of all the error areas and add them up to give the general formula:

Now we sum all the areas in the region between x = x1 and x = x2

where x2 is the upper limit and x1 is the lower limit.

Putting all these in equation (4),

 And voila!

We have our approximation:


This formula will always be more accurate than Simpson's rule and the polynomial function for a given number of iterations, since it employs both to better its approximation process.

As a brief example, we will solve the problem that we started earlier!

f(x) = sin(x) and

h(x) = 0.05621527890686895 + 0.8070522394224585.x + 0.2659311952212848.x² - 0.35406795117551315.x³ + 0.0696942723181288.x⁴ - 0.003354049885331352.x⁵


Integrating this polynomial, we get:

0.05621527890686895*x + 0.40352611971122926*x2+0.08864373174042826*x3-0.08851698779387829*x4+0.01393885446362576*x5-5.59008314221892E-4*x6

The definite integral of the polynomial beween  x = 1 and x = 3 is then easily obtained as:

h(3) - h(1) = -1.53026224468834; while its real value

Now for the summation part:

we need
 x = 1,   f(1 + 1.4)/2 = f(1.2)           h(1.2)

 x = 1.4,   f(1.4 + 1.8)/2 = f(1.6)     h(1.6)

 x = 1.8,   f(1.8 + 2.2)/2 = f(2.0)     h(2.0)

 x = 2.2,   f(2.2 + 2.6)/2 = f(2.4)     h(2.4)

 x = 2.6,   f(2.6 + 3.0)/2 = f(2.8)     h(2.8)


We present the values as coordinate pairs: (x,f,h)

(1.2, 0.932039086,  0.931961561369)


(1.6, 0.999573603,  0.999599014673)

(2.0, 0.909297427,  0.909279689992)

(2.4, 0.675463181,  0.675487289066)

(2.8, 0.334988150,  0.334918392048)

All these give our approximation as:

-1.5302948297254457

 as compared to its true value of -1.5302948024685852

Comments would be appreciated!


Thanks.

Calculating PI,

PI. That magic word that we were suddenly introduced to out of the blues in primary school.(22/7..lol)

The aim of this article is to describe some weird series I came across many years ago, while experimenting with different series.

I found that when there are certain positive integers for which:

If the integer is k: then

PI^k/constant = Summation(n ^ -k) from n=1 to n=infinity.






This series converge so rapidly..especially as the k becomes larger, that they become willing tools
in the computation of PI.



The [T, k] pairs are:


No doubt, the series gotten by substituting the first T,k values on the table is most familiar to everyone in the math. community!
Using the last pair of values on the table in the series, we get PI = 3.141592653589793 which is the highest possible accuracy available in the programming language(JAVA) which I used for the computation.


Question? how did those values come about?....not for me though..I am more of an experimental scientist than a theoretical one!



Below, I supply Java code that performs the computation:


/*
 * To change this license header, choose License Headers in Project Properties.
 * To change this template file, choose Tools | Templates
 * and open the template in the editor.
 */

package securechatengine.test;

import java.math.BigDecimal;
import java.math.MathContext;

/**
 *
 * @author JIBOYE, OLUWAGBEMIRO OLAOLUWA
 *
 */
public class TestPI {

   
    public BigDecimal calc(int terms,int power){
  BigDecimal acc = new BigDecimal(0.00,MathContext.DECIMAL128);
  BigDecimal one = new BigDecimal(1.0,MathContext.DECIMAL128);
  BigDecimal PI = new BigDecimal(Math.PI,MathContext.DECIMAL128);
        for(double i=1;i<=terms;i++){
        BigDecimal bd = new BigDecimal(i, MathContext.DECIMAL128).pow(power,MathContext.DECIMAL128);
      System.out.println( "bd = "+bd );
        acc = acc.add(one.divide(bd,MathContext.DECIMAL128));
           System.out.println( "acc = "+acc );
        }
         return PI.pow(power,MathContext.DECIMAL128).divide(acc,MathContext.DECIMAL128);   
    }
   
       public BigDecimal calcInv(int terms,int power){
  BigDecimal acc = new BigDecimal(0.00,MathContext.DECIMAL128);
  BigDecimal one = new BigDecimal(1.0,MathContext.DECIMAL128);
  BigDecimal PI = new BigDecimal(Math.PI,MathContext.DECIMAL128);
        for(double i=1;i<=terms;i++){
        BigDecimal bd = new BigDecimal(i, MathContext.DECIMAL128).pow(power,MathContext.DECIMAL128);
      System.out.println( "bd = "+bd );
        acc = i%2==0 ? acc.add(one.divide(bd,MathContext.DECIMAL128)) : acc.subtract(one.divide(bd,MathContext.DECIMAL128));
           System.out.println( "acc = "+acc );
        }
         return PI.pow(power,MathContext.DECIMAL128).divide(acc,MathContext.DECIMAL128);   
    }
   
  
   

    /**
     *
     * @param terms The number of terms of the series
     * @param power The power to use..Must be one of the magic numbers!
     * Either 2,4,6,8,10 or 14
     * @return the value of PI
     */
       public double calcPI(int terms,double power){
           double constant = getPiConstantEval(power);
         
  double acc = 0;
 

     for(double i=1;i<=terms;i++){
        double bd = Math.pow(i,power);
        acc += (1.0/bd);
        }
       double prod = constant * acc;
         return Math.pow(prod,1/power);  
    }
  
    public static void main(String[] args) {
      
         System.out.println( "Real PI = "+Math.PI );
      
       TestPI tpi = new TestPI();
     // System.out.println( tpi.calc(200, 10).toString() );
         System.out.println( "Computed PI = " + tpi.calcPI(15, 14) );
       
       
       
    }

    /**

     * This code gives the constant T for a given power

     * that we can use to compute PI.

     * You can rewrite this method to return an int or double

     * as you wish.

     */
    public static double getPiConstantEval(double power){
        switch((int)power){
          
            case 2:
                return 6.0;
            case 4:
                return 90.0;
            case 6:
                return 945.0;
            case 8:
                return 9450.0;
            case 10:
                return 93555.0;
            case 14:
                return 9121612.5;
            default:
               throw new ArithmeticException("Lol..They are magic numbers for PI dear! use them or crash it");
          
        }
    }
   
    public static BigDecimal getPiConstantEvalBigDec(double power){
        switch((int)power){
           
            case 2:
                return new BigDecimal(6, MathContext.DECIMAL128);
            case 4:
                return new BigDecimal(90, MathContext.DECIMAL128);
            case 6:
                return new BigDecimal(945, MathContext.DECIMAL128);
            case 8:
                return new BigDecimal(9450, MathContext.DECIMAL128);
            case 10:
                return new BigDecimal(93555, MathContext.DECIMAL128);
            case 14:
                return new BigDecimal(9121612.5, MathContext.DECIMAL128);
            default:
               throw new ArithmeticException("Lol..They are magic numbers for PI dear! use them or crash it");
           
        }
    }
   
}

Cache optimized Bucket sort. High speed bucket sort for floating point and integer data.

Experience is wasted if history does not repeat itself. 

 Sometime 2009, I was working on a sorting algorithm during my Student Industrial Training days and I got it working to a certain extent, but had to leave it for a while.
 It had to do with something that occurred to me then, which I called sorted clusters.

 Later, some online research showed that much work had been done in the field, and it belonged to a class of sorting algorithms called distribution sorts.

 Earlier this year, I was on a software development team and I happened to be the lead developer and we had to work on an algorithm and we agreed on sorting algorithms and then I remembered the work I did earlier, and I decided to re-work it and see what could come out of it

 That is where COBSort came from.

 It employs 3 stage cache-optimization and a formula that suggests sort indices for elements of an unsorted array to produce a high speed distribution sort.

 There are 2 flavors of this algorithm.

1. Speed Optimized COBsort

2. Memory Optimized COBsort.

 Only the first is discussed here. The difference between the two is in the extra memory needed to hold data while the sorting is been carried out.
The first needs 2 extra arrays to do this, whereas the second only needs one, making it one of the few distribution sorts to achieve this. It does this at a certain cost, however: its speed.

Its speed is usually a bit behind that of the first algorithm.  A great benefit of using only 1 extra array is that the second algorithm can handle larger input arrays than the former on limited memory systems.
 Both handle floating point data arrays and integer data arrays. Their speeds rival those of the quicksorts, mergesorts, shellsorts and the heapsorts of this world and outperform some of these algorithms in some cases.

 Java implementations of both algorithms are ready but here I give a generic algorithm for the Speed Optimized COBSort. I will post more later on.

 The time complexity of these algorithms are similar to those of the bucket sort algorithm.

 Best case = Average case = O(k*n) Worst case =O(n^2) Here is the algorithm for the Speed Optimized COBSort:

Pseudocode CobSort(inputArray,bucketSizesArray,bucketLocationsArray) 
 1. n←length[A] 
2. MAX_INSERTION_SORT←25 
3. IF n > 1000000 
4. For i = 1 to 3 do 
5. partition(array)//cache optimization //END LOOP 
6. For j = 0 to n-1 do 
7. ++bucketSizesArray[guessIndex(inputArray[j])] ///END LOOP 
8. x←0, bucketSizesAccumulator←0 
9. For k = 0 to n-1 do 
10. IF (x←bucketSizesArray[k]) > 0 
11. bucketLocationsArray[k] = bucketSizesAccumulator bucketSizesAccumulator ← bucketSizesAccumulator + x //END LOOP 
12. For i = 0 to n-1 do 
13. bucketSizesArray[bucketLocationsArray[guessIndex(inputArray[i])]++] = inputArray[i] //END LOOP 
14. For i = 1 to n-1 do 
15. inputArray[i] = bucketSizesArray[i]; 
16. For index = 0 to n-1 do 
17. start←index 
18. suggestedIndex←guessIndex(inputArray[j]) 
19. While index < n do 
20. IF suggestedIndex EQUALS guessIndex(index) 
21. index++ 
22. ELSE 
23. break //END While loop 
24. IF index GREATER THAN OR EQUALS n 
25. index = right 
26. bucketSize←index - start + 1 
27. IF bucketSize EQUALS 1 && index EQUALS n-1 
28. break 
29. ELSE IF bucketSize GREATER THAN OR. EQUALS 2 AND bucketSize LESS THAN OR EQUALS MAX_INSERTION_SORT 
30. insertionSort(inputArray,start,index) 
31. ELSE bucketSort(inputArray,start,index) //END FOR LOOP 
32. END CobSort

Thanks. Comments are welcome.

Geometry of plane areas and 2 X 2 determinants

This post gives an insight into how to compute the area of closed, regular and irregular plane areas such as rectangle, triangles, squares, trapezoidal shapes, polygons and etc. both regular and irregular.

We show a simple link between 2 X 2 determinants and plane areas, using a principle I like to call the DETERMINANT LOOP RULE. Here is a simple statement of the rule:

For a polygon having vertices (x1,y1), (x2,y2), (x3,y3),........(xn,yn), the area enclosed by the polygon is half the sum of  the 2X2 determinants of consecutive points on its vertices as we proceed in an anticlockwise loop from an initial vertex back to the same vertex.

For any 2 consecutive vertices, we form the determinant like this:

Constructing the determinant from 2 points.

To compute the area of a 3 sided plane shape, we have:

Area of a triangle with the determinant loop rule

In general, a mathematical expression of the rule is given below as:

Mathematical statement of the determinant loop rule

For example, to compute the area of the triangle whose vertices are at: (15,16), (5,5), and (20,0), we have:

Area of a triangle from first principles with the determinant loop rule

And of course, some very queer shapes can be very quickly dealt with area-wise, e.g. the very irregular shape below:

Computing the area of an irregular figure with the determinant loop rule

Pretty cool, huh?


In fact, this rule is so general that it can be applied to curves, to generate the area under a curve and so get the link between integral calculus and determinants. We could take neighboring points on a curve and let dx(the horizontal distance between them) tend to zero. We could the apply the determinant loop rule to these points and the x axis. e.g.

Integral calculus and the determinant loop rule.

As expected, the rule as stated applies directly to this diagram and computes the integral of the shape above between x1 and x2. Accuracy increase as dx tends to zero.


Further research shows that this rule is a direct link between 2X2 determinants and integral calculus.

Infact for a given dx, there is a one to one mapping between accuracies gotten here and accuracies gotten with the trapezoidal rule of integration.

It has been nice sharing knowledge  fellows. Comments would be really nice and appreciated.

Revolutionary method for preparing caustic soda and the caustic alkalis.

I came across it mid-2002 or so and since then I've been working on it, trying to understand it so that I can improve it. This technique beats every other competing technique currently in existence in terms of simplicity and is what I call the caustic soda machine of the future. Its so simple that even a youngster can set it up and have it running in a few minutes. It uses no specially designed or constructed material like asbestos or cation exchange membrane or mercury. The cell materials are cheap and readily available and the raw materials are: anolytic brine and catholytic NaOH.
In an hour, a very crude, laboratory version of the technology increased the concentration of about 50-70 ml of soda by 1.027mol/L...from 9.970 mol/L to 10.997mol/L. Despite errors in the experiment, sodium chloride impurities remained at trace levels undetected by direct electrolysis of the product. The applied voltage was about 11-12 volts and the current was between 1.2 and 1.4 amperes. The electrodes are not unusual also..I used carbon rods for both the cathode and the anode. Though I've not checked, but the current efficiency should not be too bad also. As my lab-tests come to an end on this simple, fundamental but extremely powerful invention, it is time for the invention to go to mainstream industry.

What makes the technology all the more interesting is that it may lend itself as a tool that may be useful in understanding how electric charges behave in electrolytic solutions, being a novel approach to preparing the highly valued chemical.
I wish I could go into detail about this technology here, but I cant because its not yet patented. So, till then, this is all I can reveal.

Perhaps we should start talking now.
My work has spanned calculus,mathematical tranformation of functions by rotation(in 2 and 3 dimensions),modification of Stirling's theorem to calculate factorials of all kinds of numbers,electrochemistry and .....

Hello world,this is Gbenro.The purpose of this blog is simply to tell the world about my new inventions and get people out there who are interested in them to contribute constructively.Please people of like minds,come follow this blog and let's see how far we can go.