Fourier series and square wave approximation

Fourier series is one of the most intriguing series I have met so far in mathematics. Roughly speaking it is a way to represent a periodic function using combinations of sines and cosines. There are many other fascinating topics such as the Laplace and Fourier transforms but I am new to complex analysis and techniques so I’ll go step by step!

There are many reasons why I find Fourier series fascinating but primarily, I like the fact of approximating functions using the sines, cosines and complex exponentials because of the bond between these functions.

Now, what are the ingredients we need to approximate a function with a Fourier series?

1) A periodic bounded function, with period T
2) The function should be integrable over the period

If the above conditions are met, then we can write the function as follows

Note that this is an infinite sum, the more terms you add, better the approximation. As far as the coefficients are concerned, you can calculate them using the formulas below.

For this simple test, I chose the square wave function, which could be interpreted as an on and off signal. The aim of this simple test was to check how good is the approximation. Using only 10 armonics (n=10) this is the result I got:

I would say it is pretty accurate, Python took only a few seconds to calculate it. Furthermore you can use more complex periodic functions and found similar amazing results. Here are some useful resources I used:

http://en.wikipedia.org/wiki/Fourier_series
http://en.wikipedia.org/wiki/Square_wave

If you have time, perhaps you could try plot the sawtooth wave

In the square wave Wikipedia page there are other kind of wave functions, perhaps in the future I’ll try them out too. Below is the Python code I used to plot the square wave.

 

Enjoy!

An implementation of the Tic Tac Toe game in Java

Lately I’ve been really busy, and I’ve not had that much time to spend on programming. I’ve had this implementation of the Tic Tac Toe game in my “to do list” for a very long time but sadly for one reason or another it just kept being delayed. In some sense it is easier to code than the text editor I made here, however it has not been so in practice. In particular, I can’t shake the feeling that the code of this project is redundant and sloppy. I’d be happy to be proved wrong but I guess there’s no point in arguing over such a small project.

This project, aside from getting more familiar with Java syntax and its libraries, helped me in becoming more fluent in building a GUI with Swing. There are some technical details which I’d like to discuss perhaps with some professional Java coder in case I get the chance.to talk to one (indeed networking is important).

The main Idea was to develop a Tic Tac Toe implementation and a battleship game as well, however, due to the lack of spare time I put aside the battleship game project and perhaps I will work on it in the future, although the basics of the underlying code are pretty similar to those used here (so similar that it might actually be a repetition).

Enough talk, let’s cut to the chase! Here are some screenshots, and below you can find the download link for the Java application. As usual, I’ve put into the .zip folder the source code so that you can check it out, work on it, whatever you like

How to play:
Basically you click on one of the three columns (if you click on any button on the selected column that’s also fine) and your mark (“x” for the player) is dropped into the lowest box available. Then the computer makes its move (randomly). Whoever reaches first 3 in a row.wins. The idea of the lowest box is due to the fact that I had another game in mind while developing this, it turned out as a “new” kind of Tic Tac Toe and I thought it was nice to play with. Of course by working on the source code I uploaded along with the runnable jar one could easily remove this and get back to the classic game.

Dropbox download link.

Hope this was interesting.

Taylor series with Python and Sympy

Here I am again using my beloved Python and doing maths stuff.

Today I’d like to post a short piece of code I made after a review of Taylor series I did. This script lets you input (almost) any function, provided that it can be represented using Sympy and output the Taylor series of that function up to the nth term centred at x0.

Sympy is a great module for basic symbolic mathematics, it works fine and it is really simple to use even if you are new to Python.

Here is the output of the plot function for the function sin(x) approximating up to the 9th term:

Note that in the console output the series is written backwards, however I think it could be fixed.

Here is the code I used

This code can be customized to return Taylor expansions for any function you’d like to use (of course provided it can be represented using Sympy).

Finally, the code used to generate the plot

My first serious Java project

Finally here I am with my first serious Java project: a text editor!

Yeah I know, perhaps the word “serious” is a bit of an overkill for a simple project such as a text editor. However for a Java newbie, as I consider myself to be, this was actually a middle level project. The hard part was using Swing to implement the user interface, although I find Swing to be easier to use than Tkinter (Python) once that you’ve got the basics.

So, what’s this text editor about? Basically it is a notepad (like the one that comes natively with Windows but with less options, and of course it is written in Java. Here is a screenshot:

Now, let’s talk about how to use this.

The file menu:

The file menu has 4 options: open, show current working directory, set current working directory and exit.
The open command lets you search for a .txt file to display in the user-interface.
The show and change current working directory commands are used to show and change, respectively, the directory you are working in.

You can open a file with the commands “open” and “load .txt file” or work on a new one simply by typing into the white box. When you have edited a file in the text editor, you can save it by typing a name into the bottom text field and then clicking the bottom right button. If you do not enter a name, the file will be saved with a default one. Furthermore, once you click on “save .txt data” the program asks you if you want to save the file in the current working directory. If you need to choose another location click “no” and then choose the directory to save the file in. Enter the same name of the file you are working on, to overwrite it.

The edit menu:

The edit menu basically lets you change font size (by 3 pt a click).

You can download the program (the runnable .jar file and the source code) at this address below:

https://www.dropbox.com/s/e4jhp1w9n1w6pp8/java%20text%20editor.zip?dl=0

Hope this was interesting.

Two simple Java programs

So far I’m still a newbie Java programmer and I am mostly familiar with the basics (the syntax and basic concepts) therefore I thought I could practice my little knowledge by writing two simple programs. In particular, this task helped me to get more familiar with user input, random number generation and arrays. I took the idea from a website, however I lost the bookmark therefore I can not post it.

The first program is quite simple. It asks you to guess a number within 0 and 99. Many improvements could be done to this simple program, perhaps in the future I’ll work on it.

Here it is:

The second program is a little bit more complex, it asks for an input sentence from the user and then it outputs a sentence of the same length using random part of the entered sentence.

Hope this was interesting.

Numerical integration with Python

In this short article I am going to post a simple Python script for numerical integration.

Numerical integration is a part of a family of algorithms for calculating the numerical value of a definite integral. The two simplest method for performing numerical integration are rectangle and trapezoidal rule. If you’d like to know more about these two methods you can check the wikipedia page which has been pretty helpful to me.

Now, the code below works fine if the functions you are using is not decreasing and positive, however:
-If your function is decreasing, the plot need to be adjusted to show the correct graph
-If the function is negative at some point, the code might not work, since it has not been designed to handle it. Perhaps in the future I’ll fix that.

Therefore, bottom line, square root of x is fine, logarithms (after 1),x^2,x^3(from 0) etc…

Here is an approximation with 50 points and the rectangle rule:

and the code output on the console:

Finally, below is the code I used to generate the output above:

Hope this is useful.

Tag di Technorati: maths,python,calculus,algebra,programming,integral,integrals

Volume of solids of revolution: the cone

Have you ever wondered where do all those formulas to calculate the volume of solids like a cone, a cylinder, a sphere ecc… come from? In fact they come from a simple formula and from a clever basic idea. Imagine you have a function f

Intuitively you could approximate the volume by dividing your interval into a number n of small intervals of the same size and sum up the volume of the cylinders of radius f(x) and height , then

if we use an infinite number of steps n then our sum becomes an integral:

Now that we have the intuition, we  can go on to the code. I coded an example in R and one in Matlab. The one in Matlab is shorter since it deals with the integration on its own using symbolic computation while in R I decided to code a simple approximation algorithm based on the first formula I posted above, although you could easily calculate the indefinite integral of the function foo and then use it to calculate the exact volume.

Here is the Matlab example, quick, fast and neat:

Note that f(10)=5, thus the radius of our cone is 5. Furthermore, for the sake of simplicity I’ve assumed that the tip is in 0,0. This code outputs:

Here is the R code. Of course you could have made it shorter and more neat but I like coding in R so I made something more out of it. The R program below approximates the volume using the first approach described above.

Finally, you can get a 3D cone as the one at the top of this page just by using the RGL package in R and the demo scripts.

Tag di Technorati: geometry,calculus,R,matlab,programming,volume

How to calculate the length of an arbitrary arc

Have you ever wondered if you could calculate how long is a certain curved path with a “simple formula”? Apparently you can, provided that you are good with integrals (or you can use Wolfram Alpha or Matlab). I am no mathematician and this is by no means a rigorous proof, however I think it surely delivers a first grasp on the underlying concepts.

A first approach could be use a real meter, but then you’d have to deal with the uncertainty that comes with every measurement and bla bla bla… From the deep magic forest of calculus, another approach to lines could become very handy: consider for instance an arbitrary curve S

by considering this as a limit in which the change in S approaches ds and applying the Pythagorean theorem we can write:

Now if your line is in parametric form we just divide by dt squared, if it is in the form f(x) then you can divide by dx squared and then integrate with respect to dx. For this example I am assuming S is in parametric form as such:

Dividing by dt squared yields

and after some manipulations we can integrate and get the length of our line from a to b

Suppose we’d like to calculate the length of the arc of our parametric function from 0 to 10, then assuming my calculations are correct, we’d have to calculate

Now perhaps this is easy to calculate manually, perhaps it isn’t, however it is surely fun to implement in some programming language. Sadly the sympy module in Python for symbolic maths does not work on this integral, it just seems to loop forever, perhaps that’s just my computer or I wrote something wrong. Therefore I used Matlab, which, aside from being somewhat difficult to start on my computer, works perfectly for this kind of tasks. In red is the arc whose length we’d like to calculate.

My simple Matlab script spits out as a solution the number 10.7602 which double checked on Wolfram Alpha seems to be the right solution.

To be fair my code prints out this: log(4*410^(1/2) + 81)/8 + 1640^(1/2)/4 which then computed yields the final result. Here below you can check my code

 

Tag di Technorati: maths,matlab,line,integral,programming

A plus: 2 pistons engine, get the timing right!

This is a nice add-on to my previous project: here below is a video where I simulate two pistons working together (moving a common crankshaft for instance). In order to do this effectively, one should calculate how to provide the right timing to each piston. In this example it’s quite easily done just by shifting the position of piston 2 by 180 degrees (or pi radians). A nice exercise would be trying to add more pistons.


Here is the animation on you tube:

Crankshaft connecting rod and piston mechanism simulation with Python

What is a crankshaft connecting rod and piston mechanism? It basically is a mechanical part which converts rotational motion into reciprocating motion. The applications are pretty vast, car engines are one of the most obvious I can recall right now.
It turns out the physics behind this mechanism is pretty interesting and easy to code, therefore I thought I could give it a go and try to make a simulation.
First of all we need to define the problem and solve it:



This is the basics from which we start: the crankshaft (rod ‘a’), the connecting rod (rod ‘b’) and the piston whose position is denoted with c. The motion of rod ‘a’ is a pure rotational motion while the motion of rod ‘b’ is somewhat more complex. The piston is performing a linear motion since it is constrained onto the real axis.
Now you might ask yourself why I am using the imaginary plane… It turns out that you can represent each vector (a,b and c) with a complex number, and this semplifies our problem into a more manageable system of two equations.
By using vector properties we can easily write:

The vector equation above states that the position of c is the sum of a and b as it can easily be seen on the picture. Now, we can get our real and imaginary coordinates in the plane by using Euler’s formula:

Note that every one of these parameters is a function of time and assuming that alpha(t), the length of the rods a and b are known, the system can be solved for c and beta. Furthermore by deriving the original equation we can obtain velocity and acceleration for each time t, since assuming alpha(t) is known then the derivatives of alpha(t) are known too (assuming alpha is derivable two times with respect to t).
Here are velocity and acceleration (vectors), respectively:


Now, for the sake of this example, I am assuming  however that’s not necessary. For instance one could try angular acceleration constant and so on.
Given my assumption, is known and the initial angle can be assumed later. Now we can solve our first system for and . The solution looks something like this:


Ok, so far so good, now that we have the position at each time we just need to translate everything in a language that Python can understand: I used a class, however you can easily avoid using it since it is not really necessary


Once we have coded all what is above, we can create a My_mechanism instance and call the methods, be sure not to call all the methods at once since it will not run them all. Call a method at each time:

Here below are the videos I made using the animations of the mechanism:

 

 

Hope this was interesting.

Another Excel spreadsheet: Savings with LED lights

I recently uploaded a new Excel spreadsheet where I made a calculation of how much can be saved by replacing neon tubes with LED tubes.

As you might have noticed, I do like using Excel for calculations too, and I set up a page where some of my little projects and calculations with Excel are shared. Here it is: Excel page.

Some exercises with plots and matplotlib on currencies

EDIT: Apparently, some of the prices in the .csv files I used were missing, and this caused some problems with pandas dataframe since it replaces missing values with ‘na’. Should you encounter the same problem you could check every line with an if statement or use a method to replace the na. You can check the pandas’ documentation here.

Yesterday I was really bored and had some time which could be put to good use, therefore I decided to write a quick script to plot percentage change of currencies of some countries and compare them. The code I ended up with is a bit sloppy I guess, but  that’s fine since I was primarily interested in improving my (limited, as of now) use of pandas and having fun.

First of all I gathered data from Quandl, namely the prices of the selected currencies in terms of dollars, that’s to say the value per day of every selected currency with respect to the dollar:

XYZ/USD (daily)

Gathering data from Quandl is really easy and fast using the Quandl API for Python. By the way, an API for R is available too.

I then computed the percentage change for 2 years and defined some plotting functions. Here is the main result the plotting function produces given a reference currency: it plots the percentage change for every currency (with respect to the dollar) against the percentage change of the reference currency. Some of the plotted data looks definitely weird, I wonder if I did something wrong or lost some information during the process.

Here is the code I used:

 

---

Disclaimer
This article is for educational purpose only. The author is not responsible for any consequence or loss due to inappropriate use. The article may well contain mistakes and errors. The data used might not be accurate. You should never use this article for purposes different from the educational one.

CopulaClass a Python class for using copulas: a fitting example

As I have already said in my previous post Copulalib is really user-friendly, it is difficult to write something easier, however I thought I might give it a try.

This class is built around Copulalib and since it is to be used with 2-dimensional copulas, it implements plots for data visualization and some other functions such as:
-showAvailableCopulas() a method to show visually what copulas are included in the package
-generateCopula() a method to generate the copula
-printCorrelation() this method prints out Spearman’s rho, Kendall’s tau and the fitted parameter
-getSimulatedData() this method retrieves your simulated data from the copula assuming your original data is normally distributed. It would be nice to implement some tool which could figure out the most likely distribution of your data and then use it to get the simulated observations. Perhaps in the future I’ll do it.

Furthermore, the class does not mind if you feed in python lists of numpy arrays as it turns x and y in numpy arrays. Be careful however that it does not check if your lists/arrays are of the same length.

Anyway, as for the result of the testing script below, the fitting of the Frank copula to the data seems to have been successful, our simulated data seems to fit the real quite nicely:

Originally our data was (very) approximately normally distributed, with some sort of positive correlation as you can clearly see from the plots below

Here below you can see 1000 simulated pseudo-observations from the Frank copula

Below you can find the code I used to generate this simple model:

And here are the correlation details

Copulalib: How to use copulas in Python

When dealing with copulas, R is a better option in my opinion, however, what could you do if you wish to use Python instead? There’s a good starting package called Copulalib which you can easily download here.

The package is really simple to use and very user-friendly I would say, it basically handles everything (pseudo-observations etc…) once you fed in the raw data. There is a simple example of implementation in the download page. Once the data has been fed into the function, the fitting is done automatically and the following parameters are generated:
-Spearman’s rho
-Kendall’s tau
-Theta (the parameter of the copula)

As of my understanding of the package, only Frank, Gumbel and Clayton copulas are available, and this of course could be a limitation, however it is for sure a good start. Another point which is problematic is that multidimensional copulas seem not to be supported.

Now for the real complaints: for some reason once the sample size is larger than 300 observations per variable (say 300 x and 300 y) the script raises an error saying that x and y must be of the same dimensions which is strange since they are already of the same size. Anyway maybe I did something incorrect.

Here below is the short piece of code which generated the plots of the data and of the available copulas

Next I’m going to post a class for copulas.

How to fit a copula model in R

I have been working on this topic for a great amount of time and to be honest I find R documentation not that user-friendly as the documentation for most Python modules. Anyway the fact that copulas are not the easiest model to grasp has contributed to further delays too. But mainly the lack of examples and users of these models was the biggest obstacle. Then again, I might have looked in the wrong places, if you have any good resource to suggest please feel free to leave a comment. At the bottom of this page I’ll post some links that I found very useful.

If you are new to copulas, perhaps you’d like to start with an introduction to the Gumbel copula in R here.

The package I am going to be using is the copula package, a great tool for using copulas in R. You can easily install it through R-Studio.

The dataset
For the purpose of this example I used a simple dataset of returns for stock x and y (x.txt and y.txt). You can download the dataset by clicking here. The dataset is given merely for the purpose of this example.

First of all we need to load the data and convert it into a matrix format. Optionally one can plot the data. Remember to load the copula package with library(copula)

The plot of the data

Now we have our data loaded, we can clearly see that there is some kind of positive correlation.

The next step is the fitting. In order to fit the data we need to choose a copula model. The model should be chose based on the structure of data and other factors. As a first approximation, we may say that our data shows a mild positive correlation therefore a copula which can replicate such mild correlation should be fine. Be aware that you can easily mess up with copula models and this visual approach is not always the best option. Anyway I choose to use a normal copula from the package. The fitting process anyway is identical for the other types of copula.

Let’s fit the data

Note that the data must be fed through the function pobs() which converts the real observations into pseudo observations into the unit square [0,1].
Note also that we are using the “ml” method (maximum likelihood method) however other methods are available such as “itau”.

The parameter of the fitted copula, rho, in our case is equal to 0.7387409. Let’s simulate some pseudo observations

By plotting the pseudo and simulated observations we can see how the simulation with the copula matches the pseudo observations

This particular copula might not be the best since it shows a heavy tail correlation which is not that strong in our data, however it’s a start.

Optionally at the beginning we could have plot the data with the distribution for each random variable as below

And get this beautiful representation of our original dataset

Now for the useful documentation:

Copula package official documentation:
http://cran.r-project.org/web/packages/copula/copula.pdf

R blogger article on copulas
http://www.r-bloggers.com/copulas-made-easy/

An interesting question on CrossValidated
http://stats.stackexchange.com/questions/90729/generating-values-from-copula-using-copula-package-in-r

A paper on copulas and the copula package
http://www.jstatsoft.org/v21/i04/paper

That’s all for now.