Linear programming in R

Linear programming is a technique to solve optimization problems whose constraints and outcome are represented by linear relationships.

Simply put, linear programming allows to solve problems of the following kind:

• Maximize/minimize $\hat C^T \hat X$
• Under the constraint $\hat A \hat X \leq \hat B$
• And the constraint $\hat X \geq 0$

Taylor series with Python and Sympy: Revised

More than 2 years ago I wrote a short post on Taylor series. The post featured a simple script that took a single variable function (a sine in the example), printed out the Taylor expansion up to the nth term and plotted the approximation along with the original function. As you can see on the right on the “Popular posts” bar, that post is one of the most popular and I’m told it appears among the first results on Google.

The script I wrote originally was a bit clunky, and there surely was room for improvement. Last week I received an email from a reader, Josh, who sent me an improved version of the original post.

Complex functions and flow around a cylinder

Last semester I attended a course in complex analysis and I was introduced to the basics tool used in this field. I was fascinated by almost anything that I was taught and although some topics were a little bit difficult to grasp, I am sure it was worth it all the way. Complex Analysis is great for many applications, it is a very powerful tool that enables you to solve problems in the real domain that you would have not otherwise been able to solve.

Among the many applications of complex functions, potential flow is one that took my attention, mainly because it is easygoing for a novice in this field (as I am) but also because of the nice graphs one can plot using it and their physical interpretation.

While looking up for some materials I stumbled upon this paper.

It is quite a nice piece of paper, since it gives a little bit of introduction to complex functions and how they can be used to model flow.

An ideal flow is a flow that is both incompressible and irrotational, that’s to say, if F represents the fluid velocity at each point (x,y), then the following conditions must hold in order for the flow to be ideal:

The conditions above require that the flow has no vorticity (curl) and that the volume of the fluid does not change as it flows (divergence equal to zero).

Theorem 4.5 in the paper connects the velocity vector field and complex (analytic) functions. By then integrating the velocity one can find out the streamlines of the fluid flow.

If the velocity F(z) admits a complex anti-derivative,

and

Therefore the gradient of phi is equal to F and the real part phi(x,y) defines a velocity potential for the fluid flow. The harmonic conjugate psi(x,y) to the velocity potential, in fluid mechanics, is known as the stream function.

Example (4.8). Consider the complex potential function

By making some calculations we can find out

with streamlines

asymptotically horizontal at large distances.

By using Matlab, I implemented this example and plotted some plots about our Gamma function:

 

Now we plot the velocity (gradient of the real part of Gamma) and the the streamlines (psi).

Last but not least, the contour plot of the absolute value of Gama and of its real and imaginary part.

As you can see, the stream function shows that this model can be a first raw step in the way of modelling flow around a steady cylinder. This model however is far from being perfect: among other assumptions it assumes no friction and laminar flow. Laminar flow  can be assumed when the velocity of the fluid is very low and in many practical application it is not often the case. Anyway, this simple model surely can draw some attention to complex analysis and its tools.

2D heat and wave equations on 3D graphs

While writing the scripts for the past articles I thought it might be fun to implement the 2D version of the heat and wave equations and then plot the results on a 3D graph.
As for the wave equation, Wolfram has a great page which describes the problem and explains the solution carefully describing each parameter. Here is a little animation I made using the solution they described

The code for this example is available to be downloaded here.
As far as the heat equation is concerned, I opted for a classical example where you have a hot spot in the middle of a square membrane and given an initial temperature, the heat will flow away as expected. Note that in the heat equation animation the colorbar does not change as time goes on. This can be edited by moving the colorbar into the animation function.

download the code for the heat equation.

On my youtube channel you can find many other animations of the wave equation.

The wave equation

Imagine you have an ideal string of length L and would like to find an equation that describes the oscillation of the string. Assuming the string is fixed at its ends and starts its motion in a known position f(x) the simplest assumption one can make is that the acceleration of each piece of the string is somewhat proportional to the curvature of the string as such:



We can express the considerations above in the following way:



The equation above is a partial differential equation (PDE) called the wave equation and can be used to model different phenomena such as vibrating strings and propagating waves.The constant term C has dimensions of m/s and can be interpreted as the wave speed.
It turns out that the problem above has the following general solution



The thing that strikes me about this equation is how powerful the solution is. To think about it, any function that has the argument x-ct or x+ct or a combination of both is a solution to the wave equation. This means that we can model a lot of different waves!  Furthermore, as you could probably spot, the general solution is a combination of a wave travelling to the left and one travelling to the right.
Assuming let’s try the following solution



By implementing the equation in Python for a string of length 2pi and of speed 1 m/s we obtain the following animation:

The range of animation you can run is infinite. I tried out a couple of solutions:

Shortening the string, L=pi

Click here to download the code for the above video.

Trying out f(x-ct) + g(x+ct)= cos(x-ct)**3 + cos(x+ct)**3

Click here to download the code for the above video.

The solution


yields the following:

Click here to download the code for the above video.

Note how the resulting wave (cyan) is the sum of two cosine waves travelling in opposite directions.

What happens if you input the f(x-ct) term only and set g = 0? Basically what you get is a single travelling wave. The same happens if f = 0 and g = g(x+ct).

Click here to download the code for the above video.

Hope this was entertaining and helpful!

Heat Equation part 2 a slight modification

While writing the code for the previous post I slightly modified the code in order to add 2 ‘peaks of heat’. Now, by looking at what I wrote I am not quite sure that the solution is consistent (physically speaking). At first glance we can see that there are three ‘heat spots’. If you imagine of heating the rod in two spots and cooling it in between, the situation is the following:



By setting L = 3*pi and making some other small changes we obtain:

Physically we could interpret the result as follows:
the heat flows spontaneously towards areas where the temperature is lower, therefore the temperature should go up in the middle of the rod (since it has been cooled down and has a lower temperature compared its surroundings) and go down in the sides. Assuming that the left and right boundaries of the rod are kept at constant temperature through an ideal thermostat.


Why not  checking if this solution is mathematically and physically consistent?

The Heat Equation: a Python implementation

By making some assumptions, I am going to simulate the flow of heat through an ideal rod.
Suppose you have a cylindrical rod whose ends are maintained at a fixed temperature and is heated at a certain x for a certain interval of time. Suppose that the temperature in each section with infinitesimal width dx is uniform so that we can describe the temperature in the rod using a function of only x and t.


Mathematically speaking, problem we are now facing is the following:



where k is a constant called thermal diffusivity and is different according to the different materials. By using the method of separation of variables, we can find the solution we need and by applying the initial conditions we find a particular solution for f(x) = sin(x) and L = pi
Our solution looks something like this:

Now we only need to evaluate our function at each x and t. Remember that if u(x,y) is differentiable, then:



holds. We can throw out the last term and approximate our function using the above relation since partial derivatives of u must exists and we can easily get them. Thinking about it, the second term is useless too, since we are not moving along the x axis, therefore we are left with the following:



By using matplotlib and the animation function I generated this short clip:

Here is the python implementation of the solution and the code used to graph the solution evolving with time.


Hope this was useful!

Calculate Fourier transform of a gaussian using Matlab

While wondering around in the Matlab documentation I found out there is a simple way to calculate the Fourier transform of any function using Matlab. Furthermore by using the function pretty() you can print out the result in a more friendly manner.

Let’s calculate for example the Fourier transform of the function u(x) = e^(-x^2). This calculation could also be done on paper: by applying the definition I wrote in the earlier article one finds out that:

Now, by using some properties of the transform (you can check them at http://en.wikipedia.org/wiki/Fourier_transform), we can play around and plot some functions. Here is the original function and its FT:

Let’s take the FT of e^(-(3x)^2)

It’s worth noting that when the function gets thinner, its FT gets more spread and viceversa. This phenomena is useful to understand the uncertainty principle.

This you tube video explains the link between this mathematical tools and physics.
https://www.youtube.com/watch?v=V7UNvDN_EZo

If you’d like to know more about the FT, check the wikipedia page:
http://en.wikipedia.org/wiki/Fourier_transform

Here is the Matlab code I used

Filtering an ideal signal with FFT

Hi everyone! Remember that in the last article I wrote that you can use the FFT to clean a signal from background noise? Well here is an example of signal filtering.

We start by generating a signal and then add some random noise using the random number generator in numpy. Here is the noisy signal:

Now, our signal is made up of two main frequencies: 20 and 30 Hz while the rest is mainly background noise. You can check this in the FFT of the signal:

Suppose we’d like to isolate the 20 Hz bit, that means we’d have to set to 0 each value in the Spectrum which is greater or lower than 20 Hz. Let’s block out everything outside 18 and 22 Hz and then apply the inverse FTT. Here is the signal we found (blue) compared to the exact signal we wanted to find (red). At the bottom of the graph there is the error between the two signals, notice that it is not that high, (at worst we missed 0.08, however on average we missed less).

Here is the python code I used to make this.

Have fun!

Fourier Transform: first trials

The Fourier transform is surely one of the most amazing pieces of maths applied to the sciences. For starters it is defined as follows:

The definition is fine provided that you can take the integral of u(x) over R. (Technically you should say that u(x) belongs to L1(R). But aside from these technicalities, how awesome is that definition? Calculating the FT of a function with pen and paper isn’t that straight forward at all, except in some cases. It involves mainly using complex analysis techniques. If you are familiar with calculating “standard integrals” it is rare that you can apply the same techniques here.

But how can we apply this to something? Well it turns out that if you have a signal u(x) in the time domain, its Fourier transform is the signal representation in the frequency domain. Essentially the FT enables you to decompose your time domain signal into all the frequencies that make up that signal. A first simple application is deleting noise from a signal. By checking the signal spectrum (its FT) you can choose which frequencies you need to delete from the signal in order to remove the noise. Then you just apply the inverse FT and get back your time domain signal cleaned up (more or less). Now obviously the math behind this is hairy, however this should not scare you at all, and if you are interested in these subjects I suggest you to check out some resources to get a better grasp of how all of this works.

In order to compute the FT of a signal with Python we need to use the ftt function built in into Scipy. FTT stands for Fast Fourier Transform and it is a brilliant algorithm to speed up the otherwise very long calculation of the FT of a real signal.

The signal I’ve made is a linear combination of the sine and cosine function, pretty straightforward here is its representation in the time domain:

A closer look to our signal:

Now by taking the FFT we get a list of complex numbers. Then we take the absolute value of those numbers and plot them to represent the spectrum:

Eventually, we just need to fix the lower axis and switch from bins to Hz. Since the plot is symmetrical we can take out the right half without losing information. As you can see below, we have three spikes in the frequency domain. Those three peaks corresponds to the three frequencies we set earlier in the signal (85,50,75 Hz).

 

Hope you enjoyed this, soon I’ll make a post about applying this to a real signal, such as the recording of guitar notes. Below is the code I used to plot the data above. Have fun!

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!

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

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