Friday, February 28, 2014

Plotting the time evolution of systems of differential equations of motion

Notice first that we can write each differential (dy/dx)_n as follows:
suppose y, y1, y2, y3,... are the last n values the function takes on over a certain N x values separated by a delta value.
A first derivative can be approximated like:
A second like:
And so on using divided differences of the function over a certain range.
So just picking a random relationship in x and y versus t might be:
Plugging in the approximations from above gives:

we can solve both of these for x and y, respectively, then plug in our delta for t, which I will make 1/10.:

Now we see we need to know the last two values of x (x1 and x2) and y (y1 and y2) to figure out the new x and y value.
Let's suppose the function started x1 = .5, x2 = .4, y1 = .75, y2 = 1
x = 0.5721630376
y = .7621630376
Then make these two values our x1 and y1, make the old x1 and y1: x2 and y2, and iterate. A plot of this looks like:

Sunday, February 23, 2014

Solving for a recurrence relationship of an integer sequence

Suppose you have a sequence of n numbers, it is known that there is one polynomial of degree n-1 that goes through those numbers as if they were y values when the x values were 1,2,3,n.. In the last post I discovered that those polynomials can be written as a recurrence: relationship satisfying:
-A*a(n)+B*a(n-1)-C*a(n-2)... = 0
where the coefficients A,B,C,... are row n+2 of Pascals Triangle...
But of course there are many series that if you fit a polynomial of degree n-1 through n points, that polynomial fails to accurately model the value of the n+1th point. One such series is the Fibonacci series...
I found that it's easy to find the correct recurrence formula for this series by:

Since we have eight points we look for a recurrence over the last n/2 or 4 terms, we consider that some a,b,c,d coefficients multiplying the first four terms will equal the 5th, those same a,b,c,d multiplying the 2-5th terms will equal the 6th etc. 
But looking at the recurrence formula this finds we see that d is a free variable, this means we can set it to 0 and get:
No we see that c has a variable b in it's form,  what it equals that can be set to 0 as well by making b=1.
Which gives the recurrence a(n) = 1*a(n-1)+1*a(n-2) which is the definition of the Fibonacci sequence. 
So we can look at any n numbers this way, here is one I thought of that wasn't already on the Encyclopedia of Integer sequences...
1, 3, 7, 10, 15, 23, 32, 53
If these values had been 4, -6, 4, -1 (a row on Pascals triangle with alternating signs) or it was possible through free variables to make it that, we would know that there was a cubic equation that plotted the values. Also if they could be 2, -2, 1 then a quadratic and so on. 
This one doesn't have a free variable which indicates that if there is a recurrence formula for this series it relies on at least the last 4 terms, also if there is a polynomial through the eight points it must be degree 4 or higher.  
And its predicted the next value isn't an integer, but it still gives a good guess for using 8 points. We could have also done the 9th degree polynomial through the eight points which gives:

If we had hundreds of points finding the polynomial through them gets very expensive, (exponential), while solving a set of n linear equations as in this method I've found is about O(n^3)

Another example from the Encyclopedia of Integer Sequences is the triangular numbers...
0, 1, 3, 6, 10, 15, 21, 28
Using this method this will yield the recurrence...
a(n) = 3*a(n-1)-3*a(n-2)+a(n-3)
which is the nice one found by Mark Dole in 2010. 

Wednesday, February 19, 2014

Polynomials as recurrence relations

if you have some amount of numbers you can find what the next one should be following the pattern for a degree p polynomial going through p+1 points. Like if you have -4.5, 2.3, 2, 1, you do 4*(-4.5)-6*(2.3)+4*(2)-1 and find the next number should be 11 so that a smooth curve could be drawn through the numbers if they were vertical values one apart horizontally on a graph. The 4, -6, 4, -1 are from the (p+2)th row of the Pascal's triangle.

To start consider the x^2 series starting with 0...
0, 1, 4
I noticed you can always find the next term by additions and subtractions of the initial 3 terms as follows:
4 + (4-1) + (4-1)-(1-0) = 9
This simplifies to the recurrence relationship:
A(n) = 3*(A(n-1) - A(n-2)) + A(n-3)
For example:
25 = 3*(16-9) + 4
And then I wondered what happens if you start it with other values than squares...
3, 8, 15
24 = 3*(15-8) - 3
I found the significance of this number is if you find the quadratic equation that fits these three points with any x values each a distance 1 apart, like here I randomly chose 4,5,6 as the x coordinates:

This function evaluated at the next integer 7 is 24.
The y values don't have to be positive or integers:
4.5, -3.2, 7
35.1 = 3*(7+3.2)+4.5
This function evaluated at 7 = 35.1

I also found that 3*(A(n+1)-A(n+2))+A(n+3) works equally well for finding A(n) going backwards knowing the 3 values to the right of A(n).

A similar formula can be found for cubic polynomials and it is:

 Notice the coeffiecents in both this recurrence relationship and the one for squares if the A(n) is brought to the same side read as row number P-2 where P is the power of the polynomial in Pascals triangle with alternating signs negative first, I believe this is true in general.

Now to show this one for cubes works as well:
2, -5, 7.5, 16
-3 = 4*(16)-6*(7.5)+4*(-5)-(2)

Heisenburg Uncertainty Principle makes a prime number function impossible...

First you can set up something like the following...
The way it's arranged it makes it so it is 0 when a number in the range is composite. The integer peaks are at 101, 103, 107, and 109 which are the prime numbers between 100 and 110.
So that might make you wonder if there is some function that does a better job than this one, like it has peaks of 1 at all the prime numbers and is 0 elsewhere... In the frequency domain a function like that would look like:
This is how it would look with just frequency components corresponding to the first 5 primes, you can see the peak at 1/7 in order to make a signal that repeats every 7 numbers is already having to be compressed to fit in that space. The more frequency components you would add to correspond to larger primes the narrower the frequency band has to become to fit in that space.
  This is where the Heisenburg uncertainty principle comes in, it says the narrower these bands are in the frequency domain the less certain the function is in the time domain. Meaning there couldn't be a function with information about every prime in the frequency domain that gave any information at all in the time domain. Eventually the uncertainty results in so much error as to showing which number is prime that it becomes useless.
  So there can't really be a prime number function that does something like taking on a 1 value for primes and 0 otherwise.

Outlier hypercube refinement optimization with example

I'll explain this optimization method with an example I made up, suppose you have two points and you want to find a path connecting them of length 10 when the straight line distance between them is 5. Further restrictions are that you will make your path with 6 segments, and that you want the x values of these 5 middle points to be increasing from left to right and the y values greater or equal to the starting points.
The first step of this method is to divide the search space up very roughly, To this end, I started with a further restriction that every point also be on an integer grid, and it seemed sensible to make this grid from bottom left to top right (0,0) to (5,5) in size. One can see that there would be 6*6*6*6 or 1296 possibilities to test on this grid, so I wrote a program that tested the 1296 possibilities keeping the top 5% of solutions. Meaning the 60 possibilities on this grid that came the closest in total arc length to 10 without going over 10. These 60 possibilities look like this:
This is a random selection of 13 of the top 60 solutions. Because of the constraints the x values of each solution had to be (0,1,2,3,4,5) but it was easy to see that these top 10% solutions all fell underneath the y values (0,3,4,4,3,0). So visually that cut down the search space like so:
This cuts the actually good range of the solutions to 400 possibilities or 1/3 of the original search space.
So now I looked over that range but in intervals of (1/2) instead of integer values.
Now I found that the best 10 solutions were all within 1/10000 of length 10 so I stopped. But if I wanted to be more precise the best solutions were within the red rectangles which would further cut down the search space. The black dots are one of the very best solutions found so far, there were 6 that were tied for the best:
[ ((0, 0), (1.5, 2.5), (2.5, 4.0), (3.5, 4.0), (4.5, 1.0), (5, 0)))]
which has length: 9.99856323407292

So I think this method of roughly dividing the search space, looking at the top small percentage of solutions, then refining the search space based on the extreme values of those solutions and searching that space with more refinement works well to cut down on the amount of time needed to find a very good solution.

Tuesday, February 18, 2014

Finding a zero by converging quadratic equations

I'll start with an example...
Let's say we have the quintic we want to find the zero for:
We can guess a low and high x value for example looking at the graph we might pick 1.4 and 1.5...
We find the value of the function at those x positions:
Now we can find a quadratic equation going through these above two points and a point between them at x=r whose value is 0, by solving 3 equations for the unknown coefficients. Here I'll just use Maple's interpolate command to find it...

This looks like a mess, but here we apply the idea that the closest quadratic equation to the quintic we are finding the 0 for will have to look like the quintic equation truncated to the x^2, x, and constant terms. In particular this means we can look at just the constant term of the quintic which happened to be 1 and the last term of this quadratic we just found must also have a constant term of 1. So set those equal and solve:
There might might be two roots but we want the one that is between our chosen points of 1.4 and 1.5, 1.472384325.
Now we want to evaluate the quintic at this value and see whether it is negative or positive so we know which side of the root it lies on....

Now we iterate the above steps using this new point as our new left value:
First interpolating the 3 points to find the quadratic...
Then setting the constant term of this quadratic equal to the constant term of the equation we're trying to find the root of and solving for r:
Then either accepting this as close enough to the root or continue iterating to get ever closer. 
I don't believe this is the method the computer uses but it has other ways...

Our series was 1.4, 1.472484325, 1.475958045... We were closing in fast to the approximation the computer gives of 1.476126271...

There are a lot of other ways to find roots but I thought this one was interesting. 

Sunday, February 16, 2014

Array sort revisited

Ok, made corrections to the last post. You have a list you want to sort, put it into a square array:
Sort each row individually,
Then each column individually:
Then diagonals from the top of each column toward the bottom right with wrap around.

Then the other diagonals from the top of each column towards the bottom left with wrap around...

Then those 4 steps on this matrix, 1. sort rows 2. sort columns 3. sort diagonal 4. sort reverse diagonal to get:

And I think in general it takes the base 2 log of row length number of iterations, like here it took 3 for row length 8.  
Now the list is sorted reading off by rows. 

Now this seems kind of cumbersome but the nicety of it is that each step can be carried out distributed across different processors. 
Supposing you have a number of processors equal to the square root of the number of items. This might not be as unusual as it sounds with graphics cards starting to have in the thousands of compute units.  Ordinarily n*log(n) is the best you can do on one processor, that compares to this algorithm in the sense of big(O) as:

This algorithm is the green curve. 

Parallelizable sorting algorithm (Array Sort)

**I noticed after review a counter example, and that this may only work for square matrices, the change from what I write below would be always to use a square matrix with null elements to pad it out.
There are a lot of sorting algorithms, some provably as fast on one processor as possible. But lately with the possibility of gpu and multiple core cpu architectures it is interesting to find one that can easily be split across the individual compute units acting independently. This one I call array sort which I think fits the need.
1. Start with your list of numbers or whatever is to be sorted, this is the data for my example:
In this example there are 20 items, but there can be any number. The idea for the next step is to fill an array that is as square as possible for you items.Here that would be a 4x5 array, if there had been 21 items we could use either a 4x6 array or 5x5 and fill the extra spaces with null items that would be discarded at the end. 
Now the next step is to sort each row. Since each row is sorted independently this step can be distributed across many processors for large lists. 
Now the next step is to sort by columns, this also can be distributed across many compute units. 
The third step is to sort by diagonals, where by diagonals I mean as in the linear algebra sense that each diagonal wraps around, for example the third diagonal is 3,12,1,3,15 because the third item in row one is three the fourth item in row 2 is twelve, and we wrap around to the 1 in the third row and first column and continue. This step can also be parralelized across man processors for large lists. 
Now the array can be read off by rows in sorted order. I haven't seen a case where this wasn't already sorted but if it turns out this isn't enough maybe another step where sort the top right to bottom left diagonals could be added. 
The algorithm can also be used recursively, if when we originally had a matrix that was very large, we can use the array algorithm on each row to get each of those in order etc...
This is my best attempt so far...
Suppose in our array we follow the item starting in cell [i,j] and the item starting in cell [k,l]. If they are in the same row then the larger will be put further to the right. And further sorting by column and then diagonal won't move the lesser further to the right unless it is to a row further up. If the two values don't start on the same row then they may be on the same column, in which case it will be put on a higher row and sorting by diagonals won't move the lesser to  a row further down than the greater because it would have to be greater than the number in its diagonal on the same row as the greater number but it is clearly less than that number after sorting by rows and columns.  If they start in different columns and different rows they may be on the same diagonal which will be sorted so the greater is on a further down row and further column to the right. If they aren't on the same diagonal, the lesser can't end up on a row farther down than the greater because it's been sorted by columns, and it can only move further to the right if it's on a lesser row because it's been sorted by rows. Since this is true for an arbitrary pair of numbers all the numbers will be arranged so that a number that's greater will be either on the same row and to the right of lesser numbers or on a row further down. Hence the numbers are sorted read across by rows. 
**Further idea**
Because of the recursive nature it would be sufficient if I prove by exhaustion that it works for a certain size matrix like say 4x4, then the problem could be recursively cut down until it became solving many 4x4 matrices.

Wednesday, February 12, 2014

Resonance induced chemical reaction?

As I understand NMR, when you place benzene under a very strong magnetic field you can direct different frequencies of electromagnetic radiation at the sample and there is one particular frequency where you get a resonance effect. Like in this graph:
My idea for an experiment was what if under the strong magnetic field you direct a powerful source for that particular frequency at the sample? I was thinking maybe you can put so much resonant energy in the carbon hydrogen bond that bond is stretched too far and breaks... and sometimes two carbons who have both similarly lost their hydrogen will be free to bond to each other instead of hydrogens. 
A solid sample of interconnected carbons would have great material properties, diamond is similar but might have a different crystalline structure than what you'd end up with...


I thought this was interesting...
starting with 23 and appending 3 after 3 makes a number that that is coprime to all primes less than 17.
it doesn't divide 2...
And I have a conjecture that if the number has a prime number of digits it has a good chance of being prime, but is sometimes also prime for even number of digits.

2333.... mod 3 is always 2.

it doesn't divide 5 because all multiples of 5 end in 0 or 5.
Now looking at 2333... mod 7:
23 mod 7 = 2
233 mod 7 = 2
and by induction on the division algorithm it must always be 2.

23 mod 11 = 1
233 mod 11 = 2
2333 mod 11 = 1
23333 mod 11 = 2
and also this pattern must repeat

23 mod 13 = 10
233 mod 13 = 12
2333 mod 13 = 11
23333 mod 13 = 6
233333 mod 13 = 9
2333333 mod 13 = 2
23333333 mod 13 = 10
and then repeats

mod 17 is interesting...
23333333 mod 17 is 0 so it is not prime, then it cycles through every remainder possible with 17 in a patter and
233333333333333333333333 (16 more 3's) is the next number that divides 17 evenly... and so on

Now factoring shows that:
233333 = 353*661
and 233333333 divides 29. adding 28 more 3's on the end again divides 29.

I thought maybe if there were a prime number of digits the number was prime but no. This type of number with 11, 17, and 23 and 54 digits are all prime, but 7, 13 and 19 are not. It does seem non-prime for every composite though.

Serializing an n-dimensional space?

I showed with a python program by exhaustion that:
Is unique for x, y, z an ordered triple of integers from {0..10}.
In decimals:

The closest any two values of this function over these triples is:
Which is nice because it means you could round every value to the nearest 10,000th and still have uniqueness.
I find it hard to try and figure out whether one can always map this way an ordered n-tuple {X1, X2,...Xn}
as an non-surjective injection to rational numbers with the formula {X1/P[1]+ X2/P[2]+ X3/P[N}} where P[1], ... P[N] are the N smallest primes larger than N, but I'll suggest it as a conjecture.

Sunday, February 2, 2014

Spiral Turbine

I would have thought this had been invented already if it worked but googling doesn't seem to find anything like I'm thinking...

The idea is you have a spiral extruded in the z direction off the plane inside a cylindrical housing, shaped like so...

The blue circle is a hole into the housing for an intake, and there is another hole in the housing where the arrow is. The spiral turns rapidly inside the housing (*Another idea below), and the "blade" of the spiral forces the fluid (or air) from the center to the outside along the inside of the spiral where it is compressed into a small space  between the housing and the spiral, to eventually be forced out of the hole where the arrow is.... For an incompressible fluid like water you would of course not have the housing shaped for compression.
It would take some experimenting to see what shape the spiral should be and whether it would be better to tilt the blade of the spiral off center to help guide the air but I think something along these lines would work...
I think the basic concept could be adapted for use either as actively driven for a pump, or expanding gas could turn the spiral...

*I believe it might be better if the spiral and housing are connected as one solid piece that spins within an  additional outer housing, so the effect is the same but it eliminates the friction between the spiral and the inner housing... it's hard for me to draw but the inner housing would be on a spindle attached to the top of the outer housing so it's centered and not touching the outer housing, and also have one tube for the inlet that slides in airtight solid contact with the inlet inner housing... and an output hole in the outer housing to let the pressurized fluid out...