Lagrange Interpolating Polynomial: Basic Example

The Lagrange polynomial of degree m is defined as:

The interpolation happens like

At the sampled points, the accuracy will be perfect. But at points inbetween, the accuracy of the constructed polynomial may suffer, especially for higher degree interpolations.

The runge phenomena is a manifesation of this innaccuracy, when we interpolate functions that grow with sucessive differentiation, we open ourselves up to this error, because the high order terms get large. This error appears at the boundrys, and it has been shown that using nodes with an asymptotic density to the boundrys, such as chebyshev nodes, will effectively solve this problem.

Chebyshev nodes can be constructed like:

for , the result is on an interval of

Here is a simple matlab function doing a lagrange interpolation for on to

function y2 = lagrange(x,y,x2)
for m=1:numel(x2)
	y2(m) = 0;
	for j=1:numel(x)
		l = 1;
		for k=1:numel(x)
			if k ~= j
				l = l*(x2(m)-x(k))/(x(j)-x(k));
			end

		end

	y2(m) = y2(m) + l*y(j);
	end
end

Crash Course in Perl for Data Processing: Chapter 2

Ok, now that you know some of the basics i can start to talk about useful things. This chapter is going to cover regular expressions.

Regular expressions are the most useful thing in perl, very powerful for data processing, theyre really the whole reason why im even writing this guide.

First type of regular expression: Matching

Regular expressions are enclosed in / signs.

$var =~ m/expr/; will return true if $var matches expr. Now, for some talk on how to craft “expr”

Character matches:

* – matches any character.

\w – matches whitespace character

\d – matches digits

Quantifiers:

-These specify how many times a character is going to be matched if you dont use one it will match it exactly 1 time.

. – any ammount of times

[n,m] – between n and m times

? – 0 or 1 times

Example:

$var =~ m/.*\d\d.*/;

will return true when it finds two digits somewhere in a string.

Useing matching to pull out data:

If you want to pull out characters according to a regular expression you encase the specified form in ( ) signs and you use the my() command.

my($age) = $ages =~ m/His age is: ([1,3]\d)/;

This line will fill the variable $age with the digits left in that form. We specified them as a digit between 1 and 3 characters.

Subsitutions:

if you want to subsitute all the characters of a specified type in some variable with another character you use subsitutions:

$string =~ s/a/A;

This changes all the lowercase a’s to captial A’s

Translations:

Translations are like substitutions only they work for a string, not just a single character. It allows you to substitute one string for another. Translations will only happen once per variable. so if there are multiple things in the same variable that need to be translated you are going to need to break that string up accordingly.

$string =~ tr/mike/Bob/;

if there is one mike in $string it will translate it to Bob. If there was more than one mike and you want them all translated your going to have to break the string up or hit it multiple times with the translation.

Calculating Phi

Since we are on the topic of fixed point iterative methods for calculation of irrational numbers I have decided to talk about how to calculate one of my favorite irrational numbers:

Also known as the golden ratio, is often the subject of psudoscientific conjecture and mysticism. This ratio has been used in greek architecture, under the belief that it produces aesthetically pleasing geometry. It would be used to construct golden rectangles, and as a ratio between various lengths in architecture.

is the ratio you get when you divide a line segment in two such that the ratio of the larger half to the smaller half is the same as the ratio of the whole line to the larger half. This is a pretty natural way to define it.

It also pops up in the Fibonacci sequence.

The fibonacci sequence is recursively defined as : Numbers that appear in this sequence are called the fibonacci numbers. Let denote the fibonacci number.

here is some code that will caulate it with matlab:

clear;
x(1) = 1;
x(2) = 1;
for j=3:1000
        x(j) = x(j-1) + x(j-2);
        phi = x(j)/x(j-1);
end
format long;
phi

Bojeklification for optimization and exploration of fixed point iterative methods

Bojeklification and it’s Implementation

Bojeklification is a numerical method for analysis of the convergence of iterative systems. This makes it useful for optimization of fixed point iterative methods. In order to use bojeklification one must know the attractor of the system, or in some fashion, be able to decide when a system has converged. For a fixed point this is very easy:

For some sequence that converges to a fixed point we may decide that the sequence has converged when where machine zero.

Bojeklification is when we keep track of such that the above condition is satisfied. We always have some constant in which we wish to vary, aswell as . We plot as a function of and . This gives us a simple way to see how initial conditions and some constant affects the rate of convergence.

The Fixed Point method for calculating

the following converges to when is less than about -1.5.

It is a known result, by Aimee Ross, that the optimal value is i.e. the sequence converges the fastest for this value. Also, let it be obvious that, initial values close to the attractor will also lead to fast convergence. But we wish to construct a bojeklinn of this sequence anyways, for we would wish to see if there are any other interesting things here.

Let us present the bojeklinn of this sequence:

Blue regions converge quickly, green to yellow to red converges slowly.

The thick blue line is centered at

The vertical line is at

The other line takes the form:

this is a local minimum.

Bifurcation of the Logistic Map

The logistic map is the system where bifurcation began, a man named Feigenbaund had the idea to plot the “stable values” of the system vs. values for the constant . The logistic map is recursively defined as:

The most common example of a one dimensional chaotic system, the logistic map takes on varrying properties as the value of changes. For the system converges on zero. For the attractor is the fixed point . For the attractor starts splitting into 2-cycles, 4-cycles, 8-cycles and so on into complete chaos.

A great picture illustrating this is the Feigenbaund diagram, an imaging technique otherwise known as bifurcation.

I decided to do this using the updated probabilistic method of bifurcation i developed for the lorenz attractor. Here it is:

Here is the code:

clear;
m = 1000;
M = 1000;

r = linspace(0,4,m);
N = 2000;

for k=1:numel(r)
clear x;
x(1) = 0.1;
for j=1:N
	x(j+1) = x(j)*(1-x(j))*r(k);
	X(k,j) = x(j);
end
k/(2*m)
end

draw = zeros(M,m);

xx = linspace(1,0,M);
st = round(N/10);
for j=1:m
	for k=1:M-1
		draw(k,j) = sum((xx(k) >= X(j,300:end)).*(xx(k+1) <= X(j,300:end)));
	end
(j+m)/(2*m)
end

imagesc(log10(draw+1))

Biffurcation of the Lorenz Equation

Bifurcation is a way to draw the attractor of a chaotic system. You have some iterated system with a parameter, you vary the parameter and calculate the values the system takes. The goal is to plot the “stable values”, which is most easily achievable by simply iterating the sequence for a while and plotting a number of values that the system takes after it has been allowed to stabilize.

I decided that while this method may be fine for some systems, it isn’t capable of getting a decent picture in a less predictable system with a very strange attractor. For this we want a more probabilistic approach. Some sort of probability cloud which contains the system.

How much time does the system spend on any given area? I set out to create an image that would answer this question to one extent or another. The motivation for me to develop such imaging was the lorenz equation. What is commonly referred to as the Lorenz attractor, is the lorenz equation with , and . For these values the attractor is strange. This is one of the first examples of numerical chaos.

The lorenz equation is a coupled set of differential equations:

We can varry any of these constants, , , and plot the values x, y, or z takes against them. This is a basic explanation of how standard bifurcation is done. This is also often referred to as a Feigenbaund diagram.

But to get a more complete picture, we want to take into account how often the system visits any given value.

So we divide the domain of the variable in question, x, y, or z up into discreet intervals. Then we count the number of values that fall in each bin. Just like a histogram. We tally these values up in a matrix for plotting. Every value is color coded. I use the MATLAB imagesc() command, which takes a matrix and plots it as an image, higher values on the matrix corresponding to hotter colors. I had to filter the image with a logarithm to make it more clear, without suck filtering it was very faint.

We plot the x values vs where x is the horizontal axis.

Here it is for

i apologize, the axes are not labeled correctly. Here is a zoom in on the interface between the normal attractor and advent of the strange attractor.

More interesting images to come. Here is the code i used to generate these images:

m = 2000; % rho resolution
M = 2000; % X resolution
tf = 20; %final time
R = linspace(1,60,m);

for k=1:m
clear x y z;
x(1) = 3;
y(1) = 15;
z(1) = 5;
dt = 10^-2;
pran = 10;
ray = R(k);
phy = 8/3;
N = round(tf/dt);
dt = tf/N;
for j=1:N
    x1 = x(j) + dt*(pran*(y(j)-x(j)));
    x1 = 0.75*x(j) + 0.25*(x1 + dt*(pran*(y(j)-x1)));
        x(j+1) = (x(j) + 2*(x1 + dt*(pran*(y(j)-x1))))/3;

y1 = y(j) + dt*(x(j)*(ray-z(j))-y(j));
    y1 = 0.75*y(j) + 0.25*(y1 + dt*(x(j)*(ray-z(j))-y1));
        y(j+1) = (y(j) + 2*(y1 + dt*(x(j)*(ray-z(j))-y1) ))/3;

z1 = z(j) + dt*(x(j)*y(j) - phy*z(j));
    z1 = 0.75*z(j) + 0.25*(z1 + dt*(x(j)*y(j) - phy*z1));
        z(j+1) = (z(j) + 2*(z1 + dt*(x(j)*y(j) - phy*z1)))/3;

end
k/(2*m)
Z(k,:) = z;
Y(k,:) = y;
X(k,:) = x;
end

draw = zeros(M);

xmin = min(min(X));
xmax = max(max(X));

xx = linspace(xmax,xmin,M);
st = round(N/10);
for j=1:m
        for k=1:M-1
                draw(k,j) = sum((xx(k) >= X(j,st:end)).*(xx(k+1) <= X(j,st:end))
);
        end
(j+m)/(2*m)
end
imagesc(log10(draw+1))
%save 4000x4000t.30.rho.10.2.mat draw R X Y Z
print -dpng zoomindt22.png

The Lorenz Attractor

I decided to write some code to draw trajectories in the lorenz system. Below is an image showing the trajectories of two very close initial conditions for the lorenz system.

For some time they stay close, then then diiverge

Here is the code i used

clear;
x(1) = 3;
y(1) = 15;
z(1) = 5;
x1(1) = 3.02;
y1(1) = 15.01;
z1(1) = 5.02;
dt = 10^-3;
pran = 10;
ray = 28;
phy = 8/3;
N = 20/dt;
for j=1:N
        x1(j+1) = x1(j) + dt*(pran*(y1(j)-x1(j)));
    y1(j+1) = y1(j) + dt*(x1(j)*(ray-z1(j))-y1(j));
    z1(j+1) = z1(j) + dt*(x1(j)*y1(j) - phy*z1(j));

    x(j+1) = x(j) + dt*(pran*(y(j)-x(j)));
    y(j+1) = y(j) + dt*(x(j)*(ray-z(j))-y(j));
    z(j+1) = z(j) + dt*(x(j)*y(j) - phy*z(j));
	j/N
end
plot3(x,y,z);
hold on;
plot3(x1,y1,z1,'r');

This is just a small taste of what ive been up to.

Root Finding: The Bisection Method

If there is a sign change between and and is continuous then there exists a zero for somewhere in . This is easily proved by the Intermediate value theorum which states:

If function is continuous on interval then there exists a value u between and and such that

In our case zero is between and so we know there is some value such that 

Using just this theorem it is easy to craft a numerical algorithm that will surely find the root we are looking for, given that there is only one root, and the function is continuous. We start by dividing the interval in half, and checking for a sign change on each half. By the intermediate value theorem we can conclude that the zero is in the half that has the sign change, so we discard the other half. We continue this process until our interval is of the desired accuracy, typically .

Starting with stepsize we will have a interval of size at iteration

This relationship between step and accuracy is what we mean when we talk about convergence. This algorithm crafts a sequence that converges on the value we are looking for, and the above equation tells you how fast.

Here is some code the first one is a nice clean function called like bisect(f,a,b) where f is a string containing the function and a,b is the interval. example: bisect(‘x.^2 + 3*x +1′,-1,1)

%perform the bisection method on arbitrary function
%written by Chris Bresten of Umass Dartmouth math dept.
%example: bisect(‘x.^2 + 3*x +1′,-1,1)

function z = bisect(funct,a,b);

f = inline(funct);

d = abs(b-a);
while d > 10^-15
d = abs(b-a);
c = a + ((b-a)/2);
if f(a)*f(c) < 0
b = c;
elseif f(c)*f(b) < 0
a = c;
end
if f(c) == 0
break;
end
if f(a) == 0
c = a;
break;
end
if f(b) == 0
c = b;
break;
end
end
z = c;
end

Here is a function that you call the same as the previous, this one is made as a demo to show the bisection method at work. It plots everything showing you the interval get smaller. you hit enter to make it iterate.

function z = video(funct,a,b);
f = inline(funct);
ao = a;
bo = b;
d = abs(b-a);
x = linspace(a,b);
thr = .01;
while d > 10^-15
d = abs(b-a);
c = a + ((b-a)/2);
hold off
plot(x,f(x),’k');axis([ao bo f(ao) f(bo)]);
if d < thr
x = linspace(a,b);
ao = a;
bo = b;
thr = thr/100;
end
hold on
plot(x,0,’m-’)
plot(a,0,’g.’,'MarkerSize’,50)
plot(b,0,’r.’,'MarkerSize’,50)
plot(c,0,’y*’,'MarkerSize’,50)
plot(a,f(a),’go’,'MarkerSize’,30)
plot(b,f(b),’ro’,'MarkerSize’,30)

pause
if f(a)*f(c) < 0
b = c;
elseif f(c)*f(b) < 0
a = c;
end
if f(c) == 0
break;
end
if f(a) == 0
c = a;
break;
end
if f(b) == 0
c = b;
break;
end
end
z = c;
end

Comparison of Means: Traffic at 12pm vs traffic at 8pm

For comparison of means I decided to compare the number of hits durring one hour intervals at specific times durring the day. I randomly sampled 50 days, and for each day I counted the number of requests made to the site durring the specified intervals. The intervals I chose to use were 12-1pm and 8-9pm.

The data is paired, meaning that the values in each set are connected to the corresponding values in the other set(they are collected from the same day). So to compare means it is valid to simply find the difference between the sets and compute a confidence interval for the mean of the difference.

The question we are trying to answer is weather or not there is a difference between the ammount of traffic from 12-1pm and the ammount of traffic from 8-9pm. We will compute several confidence intervals and see if any contain zero.

Here is a histogram of the differences between the ammount of traffic to the site between 12 and 1 pm and between 8 and 9 pm:

Summary Statistics are as follows:

• mean = 43.48
• median = 37.5
• skewness = 0.4640
• kurtosis = 5.0638

It appears approximately normal, low skewness, median is close to mean, kurtosis is not huge. Normal enough that im going to go ahead and do a t-interval on it.

the resulting t-interval of 95% confidence is: [21.0331 65.9269]

The confidence interval was constructed like :

where is the mean is the standard deviation and is the inverse t cdf of n degrees of freedom and confidence level

So, the confidence interval does not contain zero. I decided to try a few more samples with new intervals. Heres the rsulting t-intervals:

• [27.8191 84.9409]
• [16.2736 73.7264]
• [17.1049 66.2951]
• [51.1875 97.2925]
• [16.8320 66.8480]
• [20.2286 80.7714]

I guess we can safely conclude that there is generally more traffic from 12-1pm than from 8-9pm, perhaps people like to play on this site the most while theyre at work?

Correlagram

So I did a correlagram for the 1 hour interval counts and the one day interval counts. Check it out:

The red circles are 24, 48 and 72 hours consecutively. The peaks rest on these multiples of 24 and decay as time goes on, indicateing that the highest correlation happens when we look at it as having a 24 hour period. This can be interpreted as the level of activity on this site being periodic with a period of one day.

Now check out the 1 correlagram for 1 day intervals

You can see that the correlation decays, which really kinda shows that there is not much periodic behavior on a scale beyond 1 day. I may be able to show some periodicy by the week but that will be in another episode.

here is the definition of the correlagram for those who are interested:

where is the varience and N is half the size of the data set or less

r(k) is plotted as a function of k