Object-Oriented Programming

Python is an object-oriented programming language. The key notion is that of an object. An object consists of two things: data and functions (called methods) that work with that data. As an example, strings in Python are objects. The data of the string object is the actual characters that make up that string. The methods are things like lower, replace, and split. In Python, everything is an object. That includes not only strings and lists, but also integers, floats, and even functions themselves.

Creating your own classes.

A class is a template for objects. It contains the code for all the object’s methods.

Here is a simple example to demonstrate what a class looks like. It does not
do anything interesting.

class Esempio:
    def __init__(self, a,b):
        self.a = a
        self.b = b
    
    def add(self):
        return self.a + self.b

e = Esempio(8,6)
print(e.add())


To create a class, we use the class statement. Class names usually start with a capital.
Most classes will have a method called __init__. The underscores indicate that it is a special kind of method. It is called a constructor, and it is automatically called when someone creates a new object from your class. The constructor is usually used to set up the class’s variables. In the above program, the constructor takes two values, a and b, and assigns the class variables a and b to those values.
The first argument to every method in your class is a special variable called self. Every time your class refers to one of its variables or methods, it must precede them by self. The purpose of self is to distinguish your class’s variables and methods from other variables and functions in the program.
To create a new object from the class, you call the class name along with any values that you want to send to the constructor. You will usually want to assign it to a variable name. This is what the line e=Esempio(8,6) does.
To use the object’s methods, use the dot operator, as in e.add().

Here is a class called Analyzer that performs some simple analysis on a string. There are methods to return how many words are in the string, how many are of a given length, and how many start with a given string.

from string import punctuation

class Analizza:
    def __init__(self, s):
        for c in punctuation:
            s = s.replace(c,'')
            s = s.lower()
            self.parole = s.split()
    
    def numeri_di_parole(self):
        return len(self.parole)
    
    def inizia_con(self,s):
        return len([w for w in self.parole if w[:len(s)]==s])
    
    def numero_di_lunghezza(self, n):
        return len([w for w in self.parole if len(w)==n])


s = 'Tanto va la gatta al lardo che ci lascia lo zampino'

analizzatore = Analizza(s)

print(analizzatore.parole)
print('Numero di lettere:', analizzatore.numeri_di_parole())
print('Numero di parole che iniziano con la "c":', analizzatore.inizia_con('c'))
print('Numero di parole formato da due lettere:', analizzatore.numero_di_lunghezza(2))


A few notes about this program:

  • One reason why we would wrap this code up in a class is we can then use it a variety of different programs. It is also good just for organizing things. If all our program is doing is just analyzing some strings, then there’s not too much of a point of writing a class, but if this were to be a part of a larger program, then using a class provides a nice way to separate the Analizza code from the rest of the code. It also means that if we were to change the internals of the Analizza class, the rest of the program would not be affected as long as the interface, the way the rest of the program interacts with the class, does not change. Also, the Analizza class can be imported as-is in other programs.
  • The following line accesses a class variable: print(analizzatore.parole).

You can also change class variables. This is not always a good thing. In some cases this is convenient, but you have to be careful with it. Indiscriminate use of class variables goes against the idea of encapsulation and can lead to programming errors that are hard to fix. Some other object-oriented programming languages have a notion of public and private variables, public variables being those that anyone can access and change, and private variables being only accessible to methods within the class. In Python all variables are public, and it is up to the programmer to be responsible with them. There is a convention where you name those variables that you want to be private with a starting underscore, like _var1. This serves to let others know that this variable is internal to the class and shouldn’t be touched.

Inheritance

In object-oriented programming there is a concept called inheritance where you can create a class that builds off of another class. When you do this, the new class gets all of the variables and methods of the class it is inheriting from (called the base class). It can then define additional variables and methods that are not present in the base class, and it can also override some of the methods of the base class. That is, it can rewrite them to suit its own purposes. Here is a simple example:

class Parent:
    
    def __init__(self, a):
        self.a = a
    
    def method1(self):
        return self.a*2
    
    def method2(self):
        return self.a+'!!!'

class Child(Parent):
    def __init__(self, a, b):
        self.a = a
        self.b = b
    
    def method1(self):
        return self.a*7
    
    def method3(self):
        return self.a + self.b

p = Parent('hi')
c = Child('hi', 'bye')

print('Parent method 1: ', p.method1())
print('Parent method 2: ', p.method2())
print()
print('Child method 1: ', c.method1())
print('Child method 2: ', c.method2())
print('Child method 3: ', c.method3())


We see in the example above that the child has overridden the parent’s method1, causing it to now repeat the string seven times. The child has inherited the parent’s method2, so it can use it without having to define it. The child also adds some features to the parent class, namely a new variable b and a new method, method3.

A note about syntax: when inheriting from a class, you indicate the parent class in parentheses in the class statement.

If the child class adds some new variables, it can call the parent class’s constructor as demonstrated below. Another use is if the child class just wants to add on to one of the parent’s methods. In the example below, the child’s print_var method calls the parent’s print_var method and adds an additional line.

class Parent:
    def __init__(self, a):
        self.a = a
    
    def print_var(self):
        print("The value of this class's variables are:")
        print(self.a)

class Child(Parent):
    def __init__(self, a, b):
        Parent.__init__(self, a)
        self.b = b

def print_var(self):
    Parent.print_var(self)
    print(self.b)

Note – You can also inherit from Python built-in types, like strings (str) and lists (list), as well as any classes defined in the various modules that come with Python.
Note – Your code can inherit from more than one class at a time, though this can be a little tricky.

Uniform flow past a doublet…

Flow over a cylinder

A doublet alone does not give so much information about how it can be used to represent a practical flow pattern in aerodynamics. But let’s use our superposition powers: our doublet in a uniform flow turns out to be a very interesting flow pattern. Let’s first define a uniform horizontal flow.

u_inf = 1.0        # freestream speed

Now, we can calculate velocities and stream function values for all points in our grid. And as we now know, we can calculate them all together with one line of code per array.

u_freestream = u_inf * numpy.ones((N, N), dtype=float)
v_freestream = numpy.zeros((N, N), dtype=float)

psi_freestream = u_inf * Y

Below, the stream function of the flow created by superposition of a doublet in a free stream is obtained by simple addition.

The plot shows that this pattern can represent the flow around a cylinder with center at the location of the doublet. All the streamlines remaining outside the cylinder originated from the uniform flow. All the streamlines inside the cylinder can be ignored and this area assumed to be a solid object. This will turn out to be more useful than you may think.

# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')

# calculate the stagnation points
x_stagn1, y_stagn1 = +math.sqrt(kappa / (2 * math.pi * u_inf)), 0.0
x_stagn2, y_stagn2 = -math.sqrt(kappa / (2 * math.pi * u_inf)), 0.0

# display the stagnation points
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');

Bernoulli’s equation and the pressure coefficient

A very useful measurement of a flow around a body is the coefficient of pressureCp. To evaluate the pressure coefficient, we apply Bernoulli’s equation for ideal flow and with simple mathematical steps:

\(C_{p}=1-\left(\frac{U}{U_{\infty}}\right)^{2}\)

In an incompressible flow, Cp=1 at a stagnation point. Let’s plot the pressure coefficient in the whole domain.

# compute the pressure coefficient field
cp = 1.0 - (u**2 + v**2) / u_inf**2

# plot the pressure coefficient field
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(1.1 * width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
contf = pyplot.contourf(X, Y, cp,
                        levels=numpy.linspace(-2.0, 1.0, 100), extend='both')
cbar = pyplot.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyplot.scatter(x_doublet, y_doublet,
               color='#CD2305', s=80, marker='o')
pyplot.contour(X,Y,psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyplot.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');
pyplot.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')

Below, we report the complete Python code:

import math as mt
import numpy as np
from matplotlib import pyplot as pyp

N = 50                                # Number of points in each direction
x_start, x_end = -2.0, 2.0            # x-direction boundaries
y_start, y_end = -1.0, 1.0            # y-direction boundaries
x = np.linspace(x_start, x_end, N)    # creates a 1D-array for x
y = np.linspace(y_start, y_end, N)    # creates a 1D-array for y
X, Y = np.meshgrid(x, y)              # generates a mesh grid

kappa = 1.0                        # strength of the doublet
x_doublet, y_doublet = 0.0, 0.0    # location of the doublet

def get_velocity_doublet(strength, xd, yd, X, Y):
    
    u = (- strength / (2 * mt.pi) *
         ((X - xd)**2 - (Y - yd)**2) /
         ((X - xd)**2 + (Y - yd)**2)**2)
    v = (- strength / (2 * mt.pi) *
         2 * (X - xd) * (Y - yd) /
         ((X - xd)**2 + (Y - yd)**2)**2)
    
    return u, v

def get_stream_function_doublet(strength, xd, yd, X, Y):
    
    psi = - strength / (2 * mt.pi) * (Y - yd) / ((X - xd)**2 + (Y - yd)**2)
    
    return psi

# compute the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)

# compute the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
pyp.streamplot(X, Y, u_doublet, v_doublet,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyp.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o');



u_inf = 1.0        # freestream speed

u_freestream = u_inf * np.ones((N, N), dtype=float)
v_freestream = np.zeros((N, N), dtype=float)

psi_freestream = u_inf * Y

# superposition of the doublet on the freestream flow
u = u_freestream + u_doublet
v = v_freestream + v_doublet
psi = psi_freestream + psi_doublet

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
pyp.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyp.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyp.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o')

# calculate the stagnation points
x_stagn1, y_stagn1 = +mt.sqrt(kappa / (2 * mt.pi * u_inf)), 0.0
x_stagn2, y_stagn2 = -mt.sqrt(kappa / (2 * mt.pi * u_inf)), 0.0

# display the stagnation points
pyp.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');

# compute the pressure coefficient field
cp = 1.0 - (u**2 + v**2) / u_inf**2

# plot the pressure coefficient field
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyp.figure(figsize=(1.1 * width, height))
pyp.xlabel('x', fontsize=16)
pyp.ylabel('y', fontsize=16)
pyp.xlim(x_start, x_end)
pyp.ylim(y_start, y_end)
contf = pyp.contourf(X, Y, cp,
                        levels=np.linspace(-2.0, 1.0, 100), extend='both')
cbar = pyp.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyp.scatter(x_doublet, y_doublet,
               color='#CD2305', s=80, marker='o')
pyp.contour(X,Y,psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid')
pyp.scatter([x_stagn1, x_stagn2], [y_stagn1, y_stagn2],
               color='g', s=80, marker='o');
pyp.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')


pyp.show()

Solving 2D Heat Equation…

Before we do the Python code, let’s talk about the heat equation and finite-difference method. Heat equation is basically a partial differential equation, it is

If we want to solve it in 2D (Cartesian), we can write the heat equation above like this:

where u is the quantity that we want to know, t is for temporal variable, x and y are for spatial variables, and α is diffusivity constant. So basically we want to find the solution u everywhere in x and y, and over time t.

We can write the heat equation above using finite-difference method like this:

If we arrange the equation above by taking Δx = Δy, we get this final equation:

where

\(\gamma = \alpha \frac{{\Delta t}}{{\Delta {x^2}}}\)

We use explicit method to get the solution for the heat equation, so it will be numerically stable whenever

\(\Delta t \leq \frac{{\Delta {x^2}}}{{4\alpha }}\)

Everything is ready. Now we can solve the original heat equation approximated by algebraic equation above, which is computer-friendly.

Let’s suppose a thin square plate with the side of 50 unit length. The temperature everywhere inside the plate is originally 0 degree (at t = 0), let’s see the diagram below:

For our model, let’s take Δx = 1 and α = 2.0.
Now we can use Python code to solve this problem numerically to see the temperature everywhere (denoted by i and j) and over time (denoted by k). Let’s first import all of the necessary libraries, and then set up the boundary and initial conditions.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import FuncAnimation

print("2D heat equation solver")

plate_length = 50
max_iter_time = 750

alpha = 2
delta_x = 1

delta_t = (delta_x ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_x ** 2)

# Initialize solution: the grid of u(k, i, j)
u = np.empty((max_iter_time, plate_length, plate_length))

# Initial condition everywhere inside the grid
u_initial = 0

# Boundary conditions
u_top = 100.0
u_left = 0.0
u_bottom = 0.0
u_right = 0.0

# Set the initial condition
u.fill(u_initial)

# Set the boundary conditions
u[:, (plate_length-1):, :] = u_top
u[:, :, :1] = u_left
u[:, :1, 1:] = u_bottom
u[:, :, (plate_length-1):] = u_right

We’ve set up the initial and boundary conditions, let’s write the calculation function based on finite-difference method that we’ve derived above.

def calculate(u):
    for k in range(0, max_iter_time-1, 1):
        for i in range(1, plate_length-1, delta_x):
            for j in range(1, plate_length-1, delta_x):
                u[k + 1, i, j] = gamma * (u[k][i+1][j] + u[k][i-1][j] + u[k][i][j+1] + u[k][i][j-1] - 4*u[k][i][j]) + u[k][i][j]

    return u

Let’s prepare the plot function so we can visualize the solution (for each k) as a heat map. We use Matplotlib library, it’s easy to use.

def plotheatmap(u_k, k):
    # Clear the current plot figure
    plt.clf()

    plt.title(f"Temperature at t = {k*delta_t:.3f} unit time")
    plt.xlabel("x")
    plt.ylabel("y")

    # This is to plot u_k (u at time-step k)
    plt.pcolormesh(u_k, cmap=plt.cm.jet, vmin=0, vmax=100)
    plt.colorbar()

    return plt

One more thing that we need is to animate the result because we want to see the temperature points inside the plate change over time. So let’s create the function to animate the solution.

def animate(k):
  plotheatmap(u[k], k)

anim = animation.FuncAnimation(plt.figure(), animate, interval=1, frames=max_iter_time, repeat=False)
anim.save("heat_equation_solution.gif")

That’s it! And here’s the result:

…a doublet.

Think about the Source & Sink again, and now imagine that you are looking at this flow pattern from very far away. The streamlines that are between the source and the sink will be very short, from this vantage point. And the other streamlines will start looking like two groups of circles, tangent at the origin. If you look from far enough away, the distance between source and sink approaches zero, and the pattern you see is called a doublet.

Let’s see what this looks like. First, load our favorite libraries.

First, consider a source of strength \(\Lambda\) at (−l2,0) and a sink of opposite strength located at (l2,0). At any point P in the flow, the stream funcion is

\(\psi(x, y)=\frac{\Lambda}{2 \pi}\left(\theta_{1}-\theta_{2}\right)=-\frac{\Lambda}{2 \pi} \Delta \theta \)

where \(\Delta \theta = \theta_2 – \theta_1\). Let the distance l between the two singularities approach zero while the strength magnitude is increasing so that the product \(\Lambda l\) remains constant. In the limit, this flow pattern is a doublet and we define its strength by \(k=l \Lambda\). The stream function for a doublet is obtained from previous equation as follows:

\(\psi \left( {x,y} \right) = \mathop {\lim }\limits_{l \to 0} \left( { – \frac{\Lambda }{{2\pi }}d\theta } \right){\text{ and }}k = l\Lambda = {\text{cost}}\)

where in the limit \(\Delta \theta \to d\theta \to 0\).

Considering the case where is infinitesimal, we deduce from the figure above that:

Hence the stream function becomes:

i.e.

\(\psi(r, \theta)=-\frac{\kappa}{2 \pi} \frac{\sin \theta}{r}\)

In Cartesian coordinates, a doublet located at the origin has the stream function

\(\psi(x, y)=-\frac{\kappa}{2 \pi} \frac{y}{x^{2}+y^{2}}\)

from which we can derive the velocity components

Now we have done the math, it is time to code and visualize what the streamlines look like. We start by creating a mesh grid.

N = 50                                # Number of points in each direction
x_start, x_end = -2.0, 2.0            # x-direction boundaries
y_start, y_end = -1.0, 1.0            # y-direction boundaries
x = numpy.linspace(x_start, x_end, N)    # creates a 1D-array for x
y = numpy.linspace(y_start, y_end, N)    # creates a 1D-array for y
X, Y = numpy.meshgrid(x, y)              # generates a mesh grid

We consider a doublet of strength κ=1.0 located at the origin.

kappa = 1.0                        # strength of the doublet
x_doublet, y_doublet = 0.0, 0.0    # location of the doublet

We play smart by defining functions to calculate the stream function and the velocity components that could be re-used if we decide to insert more than one doublet in our domain.

Once the functions have been defined, we call them using the parameters of the doublet: its strength kappa and its location x_doublet, y_doublet.

# compute the velocity field on the mesh grid
u_doublet, v_doublet = get_velocity_doublet(kappa, x_doublet, y_doublet, X, Y)

# compute the stream-function on the mesh grid
psi_doublet = get_stream_function_doublet(kappa, x_doublet, y_doublet, X, Y)

We are ready to do a nice visualization.

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u_doublet, v_doublet,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(x_doublet, y_doublet, color='#CD2305', s=80, marker='o');

pyplot.show()

…in a Freestream

Starting from what is written here, we will try to modularize our code using function definitions. This will make things easier to manage.

We start by importing the libraries that we will be using in our code: the NumPy array library, the Matplotlib plotting library and the mathematical functions in the math module.

import numpy as np
import math as mt
from matplotlib import pyplot
from matplotlib import cm

To visualize the streamlines, we need to create a grid of points where we’ll compute the velocity.

We have our set of points now, and the two arrays X and Y contain their x– and y– coordinates (respectively) of every point on the rectangular grid.

Source in a uniform flow

We will first superimpose a source on a uniform flow and see what happens.

The streamlines of a freestream with speed \(U_\infty\) and angle of attack α are given by:

\(\psi_{\text {freestream }}(x, y)=U_{\infty}(y \cos \alpha-x \sin \alpha)\)

Note: the streamlines are all straight, parallel lines that make an angle α with the x-axis. If the flow is completely horizontal, ψ=Uy. Integrate, and you get that u=U∞ and v=0.

Let’s write some code that will fill the arrays containing the u-velocity, the v-velocity and the stream function of a uniform horizontal flow (U,α=0), on every point of our grid. Note the handy NumPy functions ones(), which creates a new array and fills it up with the value 1 everywhere, and zeros(), which creates an array filled with 0.

The stream function of a source flow located at (xsource,ysource) is:

\(\psi_{\text {source }}(x, y)=\frac{\sigma}{2 \pi} \arctan \left(\frac{y-y_{\text {source }}}{x-x_{\text {source }}}\right)\)

and the velocity components are:

And remember that the stream function and velocity field of a source and a sink are exactly the same except one has positive strength while the other has negative strength.

We can write functions that serve a double purpose: with σ positive, they give the velocity and stream function of a source; with σ negative, they give them for a sink.

def get_velocity(strength, xs, ys, X, Y):

    u = strength / (2 * np.pi) * (X - xs) / ((X - xs)**2 + (Y - ys)**2)
    v = strength / (2 * np.pi) * (Y - ys) / ((X - xs)**2 + (Y - ys)**2)
    
    return u, v

Note that the output of the function consists of two arrays: u and v. They are calculated inside the function, which is indicated by the indentation of the lines after the colon. The final line indicates with the return keyword that the arrays u, v are sent back to the statement that called the function.

Similarly, we define another function to compute the stream-function of the singularity (source or sink) on the mesh grid, and call it get_stream_function().

def get_stream_function(strength, xs, ys, X, Y):
 
    psi = strength / (2 * np.pi) * np.arctan2((Y - ys), (X - xs))
    
    return psi

Let’s use our brand new functions for the first time:

strength_source = 5.0            # strength of the source
x_source, y_source = -1.0, 0.0   # location of the source

# compute the velocity field
u_source, v_source = get_velocity(strength_source, x_source, y_source, X, Y)

# compute the stream-function
psi_source = get_stream_function(strength_source, x_source, y_source, X, Y)

Let’s again use our superposition powers. The streamlines of the combination of a freestream and a source flow are:

\(\psi=\psi_{\text {freestream }}+\psi_{\text {source }}=U_{\infty} y+\frac{\sigma}{2 \pi} \arctan \left(\frac{y-y_{\text {source }}}{x-x_{\text {source }}}\right)\)

And since differentiation is linear, the velocity field induced by the new flow pattern is simply the sum of the freestream velocity field and the source velocity field:

The stagnation points in the flow are points where the velocity is zero. To find their location, we solve the following equations:

\(u=0, v=0\)

which leads to:

The streamline containing the stagnation point is called the dividing streamline. It separates the fluid coming from the freestream and the fluid radiating from the source flow. On the streamline plot, we’ll add a red curve to show the dividing streamline, and we’ll use the contour() function for that.

We will also draw a red circle to show the location of the stagnation point, using the scatter() function.

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.grid(True)
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v, density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter(x_source, y_source, color='#CD2305', s=80, marker='o')

# calculate the stagnation point
x_stagnation = x_source - strength_source / (2 * np.pi * u_inf)
y_stagnation = y_source

# display the stagnation point
pyplot.scatter(x_stagnation, y_stagnation, color='g', s=80, marker='o')

# display the dividing streamline
pyplot.contour(X, Y, psi, 
               levels=[-strength_source / 2, strength_source / 2], 
               colors='#CD2305', linewidths=2, linestyles='solid');

pyplot.show()

If we ignore the flow inside the dividing streamline, we can consider that a solid body. In fact, this body has a name: it is called a Rankine half body.

Source-sink pair in a uniform flow

Now we can add a sink to our flow pattern without too much extra coding.

strength_sink = -5.0        # strength of the sink
x_sink, y_sink = 1.0, 0.0   # location of the sink

# compute the velocity field on the mesh grid
u_sink, v_sink = get_velocity(strength_sink, x_sink, y_sink, X, Y)

# compute the stream-function on the grid mesh
psi_sink = get_stream_function(strength_sink, x_sink, y_sink, X, Y)

The superposition of the freestream, the source and the sink is just a simple addition.

# superposition of a source and a sink on the freestream
u = u_freestream + u_source + u_sink
v = v_freestream + v_source + v_sink
psi = psi_freestream + psi_source + psi_sink

# plot the streamlines
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
pyplot.streamplot(X, Y, u, v,
                  density=2, linewidth=1, arrowsize=1, arrowstyle='->')
pyplot.scatter([x_source, x_sink], [y_source, y_sink],
               color='#CD2305', s=80, marker='o')
pyplot.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid');

We can can look at the elliptical closed streamline as a solid surface and imagine that this is the flow around an egg-shaped object. It is called a Rankine oval.

Bernoulli’s equation and the pressure coefficient

A very useful measurement of a flow around a body is the coefficient of pressure Cp. To evaluate the pressure coefficient, we apply Bernoulli’s equation for an incompressible flow:

We define the pressure coefficient in the following way:

i.e.,

Note that in an incompressible flow, Cp=1 at a stagnation point.

# compute the pressure coefficient field
cp = 1.0 - (u**2 + v**2) / u_inf**2

# plot the pressure coefficient field
width = 10
height = (y_end - y_start) / (x_end - x_start) * width
pyplot.figure(figsize=(1.1 * width, height))
pyplot.xlabel('x', fontsize=16)
pyplot.ylabel('y', fontsize=16)
pyplot.xlim(x_start, x_end)
pyplot.ylim(y_start, y_end)
contf = pyplot.contourf(X, Y, cp,
                        levels=numpy.linspace(-2.0, 1.0, 100), extend='both', cmap=cm.gray)
cbar = pyplot.colorbar(contf)
cbar.set_label('$C_p$', fontsize=16)
cbar.set_ticks([-2.0, -1.0, 0.0, 1.0])
pyplot.scatter([x_source, x_sink], [y_source, y_sink],
               color='#CD2305', s=80, marker='o')
pyplot.contour(X, Y, psi,
               levels=[0.], colors='#CD2305', linewidths=2, linestyles='solid');

Equazioni differenziali ordinarie…

L’idea è di scrivere un codice per risolvere l’equazione differenziale dell’oscillatore armonico generica, implementando i metodi di : Eulero, Eulero-Cromer e Verlet.

Gli schemi di integrazione sono:

  • Eulero: \(v_{n+1}=v_{n}+\tau F_{n} ; x_{n+1}=x_{n}+\tau v_{n}\)
  • Eulero-Cromer: \(v_{n+1}=v_{n}+\tau F_{n} ; x_{n+1}=x_{n}+\tau v_{n+1}\)
  • Verlet: \(x_{n+1}=x_{n}+\tau v_{n}+F_{n} \frac{\tau^{2}}{2} ; v_{n+1}=v_{n}+\frac{\tau^{2}}{2}\left(F_{n}+F_{n+1}\right) \)

dove con \(\tau\) indichiamo il passo temporale e con \(F_n\) indichiamo l’accelerazione, o la forza diviso la massa, a cui è sottoposto il nostro sistema. Nel nostro caso avremo:

\(\frac{d^{2} \theta}{d t^{2}}=-\frac{g}{l} \sin \theta=-\omega^{2} \sin \theta\)

La traduzione in codice di questi schemi è immediata. Vediamola di seguito e poi ci costruiamo il programma intorno:

#Accelerazione
F0 = -g/l*sin(x0)

#Eulero
v1 = v0+tau*F0
x1 = x0+tau*v0

#Eulero-Cromer
v1 = v0+tau*F0
x1 = x0+tau*v1

#velocity-Verlet
x1 = x0+v0*tau+0.5*F0*tau**2
F1 = -g/l*sin(x1)
v1 = v0+0.5*tau*(F0+F1)

Nel codice, facciamo uso di due variabili una denominata con 0 l’altra con 1 che si riferiscono al tempo \(n\) e \(n+1\) rispettivamente. Se ora si definiscono le costanti necessarie (g, l, τ) e le condizioni iniziali (x0,v0), la serie di righe scritte sopra ci permettono di calcolare la posizione e la velocità dopo un tempo τ, con il metodo che desideriamo.
Vediamo di seguito tali definizioni:

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi

g = 1
l = 1
x0 = np.deg2rad(5)
v0 = 0.0
tau = 0.01

Facciamo notare l’uso della funzione deg2rad(alpha) che permette di convertire l’angolo alpha espresso in gradi sessagesimali in radianti. Questo è necessario perché la funzione sin(alpha) che troviamo per calcolare la forza prende in input un angolo alpha in radianti. Esiste anche la funzione inversa rad2deg().

Ovviamente vorremo prolungare la nostra simulazione, calcolando la posizione e la velocità del nostro oscillatore, per un periodo di tempo più lungo di in solo passo τ, per esempio per T = N · τ, con N a piacere. E facile capire che dovremo ripetere lo schema di integrazione scelto per un numero N di volte, mediante un ciclo for,
come vediamo di seguito:

N = 3000
for t_step in range(N):
    #Accelerazione
    F0 = -g/l*sin(x0)
    #Eulero-Cromer
    v1 = v0+tau*F0
    x1 = x0+tau*v1
    x0 = x1
    v0 = v1

Nelle ultime due righe ci ricordiamo di aggiornare la variabile della posizione e della velocità denominata con 0 con quelle appena calcolate, denominate con 1.
A questo punto il nostro programma calcola la posizione e la velocità dell’oscillatore per gli N istanti di durata τ, ma è ancora poco utile, perché dopo l’uscita dal ciclo for avremo a disposizione solo l’ultimo valori di posizione e velocità. Sarebbe utile quindi registrare l’intera dinamica per poi studiarne l’andamento mediante un grafico, per esempio. Facciamo in questo modo:

• Nelle righe iniziali del programma, quindi prima del ciclo for, creiamo due array:

N = 3000

x_vet = np.zeros(N+1)
v_vet = np.zeros(N+1)

• registriamo nel primo elemento degli array le condizioni iniziali:

x_vet[0] = x0
v_vet[0] = v0

• all’interno del ciclo for registriamo ogni nuovo valore di x e v all’interno dei rispettivi array:

for timestep in range(N):
    #Accelerazione
    F0 = -g/l*sin(x0)

    x1 = x0+v0*tau+0.5*F0*tau**2
    F1 = -g/l*sin(x1)
    
    v1 = v0+0.5*tau*(F0+F1)
    
    x0 = x1
    v0 = v1
    x_vet[timestep+1] = x0
    v_vet[timestep+1] = v0

Ora una volta eseguite le righe di codice precedenti potremo fare il grafico della posizione (velocità) in funzione del tempo, creando prima un array con tutti gli istanti temporali usati nella integrazione, mediante la funzione linspace():

time_vet = np.linspace (0,N*tau, N+1)

plt.xlabel('Tempo (sec)')
plt.ylabel('Ampiezza (rad)')

plt.plot(time_vet, x_vet)
plt.show()

Ecco il codice completo:

import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi


g = 1
l = 1
x0 = np.deg2rad(5)
v0 = 0.0
tau = 0.01


N = 3000

x_vet = np.zeros(N+1)
v_vet = np.zeros(N+1)

x_vet[0] = x0
v_vet[0] = v0

for timestep in range(N):
    #Accelerazione
    F0 = -g/l*sin(x0)

    x1 = x0+v0*tau+0.5*F0*tau**2
    F1 = -g/l*sin(x1)
    
    v1 = v0+0.5*tau*(F0+F1)
    
    x0 = x1
    v0 = v1
    x_vet[timestep+1] = x0
    v_vet[timestep+1] = v0

time_vet = np.linspace (0,N*tau, N+1)

plt.xlabel('Tempo (sec)')
plt.ylabel('Ampiezza (rad)')

plt.plot(time_vet, x_vet)
plt.show()

e questo è il diagramma ottenuto:

Mathematical induction…

There are many techniques to prove various statements. The simplest ones are direct proof, where you show that if the first statement (P) is true then the second statement (Q) is also true; the proof by contrapositive, where you assume that the negation of the second statement (~Q) is true and show that this forces the negation of the first statement (~P) to be true as well. Mathematical induction is a technique used to prove that an all statements of an infinite list of statements (instead of just two like in the first two techniques mentioned) are all true.

To motivate the discussion, let’s first examine the kinds of statements that induction is used to prove.
Consider the following statement.

Conjecture. The sum of the first \(n\) odd natural numbers equals \(n^2\).

The following table illustrates what this conjecture says. Each row is headed by a natural number \(n\), followed by the sum of the first \(n\) odd
natural numbers, followed by \(n^2\).

Note that in the first five lines of the table, the sum of the first \(n\) odd numbers really does add up to \(n^2\). Notice also that these first five lines indicate that the \(n\)th odd natural number (the last number in each sum)
is \(2n-1\). (For instance, when \(n=2\), the second odd natural number is \(2 \cdot 2−1 =3\); when \(n=3\), the third odd natural number is \(2 \cdot 3−1 =5\), etc.)
The table raises a question. Does the sum \(1+3+5+7+\cdots+(2n−1)\) really always equal \(n^2\)? In other words, is the conjecture true?
Let’s rephrase this as follows. For each natural number \(n\) (i.e., for each
line of the table), we have a statement \(S_n\), as follows:

Our question is: Are all of these statements true?

Mathematical induction is designed to answer just this kind of question.
It is used when we have a set of statements \(S_1,S_2,S_3,\cdots ,S_n,\cdots \), and weneed to prove that they are all true. The method is really quite simple.

In this setup, the first step (1) is called the basis step. Because \(S1\) is usually a very simple statement, the basis step is often quite easy to do.
The second step (2) is called the inductive step. In the inductive step direct proof is most often used to prove \(S_k \Rightarrow S_{k+1}\), so this step is usually carried out by assuming \(S_k\) is true and showing this forces \(S_{k+1}\) to be true.
The assumption that \(S_k\) is true is called the inductive hypothesis.

Now let’s apply this technique to our original conjecture that the sum of the first \(n\) odd natural numbers equals \(n^2\). Our goal is to show that for each \(n \in \mathbf{N}\), the statement \(S_n : 1+3+5+7+\cdots+(2n−1) = n^2\) is true. Before getting started, observe that \(S_k\) is obtained from \(S_n\) by plugging \(k\) in for \(n\).
Thus \(S_k\) is the statement \(S_k : 1+3+5+7+\cdots+(2k −1)=k^2\). Also, we get \(S_{k+1}\) by plugging in \(k +1\) for \(n\), so that \(S_{k+1} : 1+3+5+7+\cdots+(2(k +1)−1) =(k +1)^2\).

More:

Using AnimatPlot for Animating Graphs & Plots

Data Visualization helps in understanding different patterns, associations, visual insights from the data, etc. It is important because it uncovers the mystery behind the data tables in form of charts, graphs, and plots. There are tons of python libraries that help in visualizing the data like Matplotlib, Seaborn, etc.

AnimatPlot is an open-source python library that is built on top of Matplotlib and is used for creating highly interactive animated plots. Now, we will explore some of the functionalities that AnimatPlot provides.

Installing required libraries

We will start by installing AnimatPlot using pip. The command given below will do that.

pip install animatplot

Animatplot is built on the concept of blocks. We’ll start by animating a Line block.

First we need some imports.

Interactivity is not available in the static docs. Run the code locally to get interactivity.

So, we will use JupyterLab for our simulations.

We will import the required libraries for creating animated plots.

We will animate the function:

\(y=\cos (2 \pi(x+t)) \text { over the range } x=[0,1] \text {, and } t=[0,1]\)

Let’s generate the data:

In order to tell animatplot how to animate the data, we must pass it into a block. By default, the Line block will consider each of the rows in a 2D array to be a line at a different point in time.

We then pass a list of all our blocks into an Animation, and show the animation.

We’ll use the same data to make a new animation with interactive controls.

block = amp.blocks.Line(X, Y)
anim = amp.Animation([block])

anim.controls() # creates a timeline_slider and a play/pause toggle
anim.save_gif('images/line2') # save animation for docs
plt.show()

Displaying the Time

The above animation didn’t display the time properly because we didn’t tell animatplot what the values of time are. Instead it displayed the frame number. We can simply pass our values of time into our call to Animation.

block = amp.blocks.Line(X, Y)
anim = amp.Animation([block], t) # pass in the time values

anim.controls()
anim.save_gif('line3') # save animation for docs
plt.show()

Similarly, now we will add more data and create multiple plots in a single chart.

and once again:

Insieme di Mandelbrot…

L’insieme di Mandelbrot o frattale di Mandelbrot è l’insieme dei numeri complessi \(c \in \displaystyle{C}\) per i quali la successione definita da:

{\displaystyle {\begin{cases}z_{0}=0\\z_{n+1}=z_{n}^{2}+c\end{cases}}}

è limitata.

Nonostante la semplicità della definizione, l’insieme ha una forma complessa il cui contorno è un frattale.

Qui, la funzione per la costruzione della matrice da diagrammare:

def plotter(n, vallim, itlim, xi,xs,yi,ys):
    image = np.full((n,n),0)
    Xp = np.linspace(xi,xs,n)
    Yp = np.linspace(yi,ys,n)
    for x in range(n):
        for y in range(n):
            image[y][x] = divergetest(complex(Xp[x],Yp[y]),vallim, itlim)
    return image

Di seguito, invece, la routine per verificare la convergenza della successione:

def divergetest(c, vallim, itlim):
    z = c
    i = 0
    while i < itlim and np.sqrt((z.real)**2 + (z.imag)**2)< vallim:
        z = z**2 + c
        i = i + 1
    return i

Il programma principale:

n = 5000
xi,xs,yi,ys = -3,1,-2,2
vallim = 5
itlim = 30
mandelbrot = plotter(n,vallim,itlim,xi,xs,yi,ys)

plt.figure(1,figsize=[8,8])
plt.title('Mandelbrot Set')
plt.imshow(mandelbrot,extent=[-3,1,-2,2])
plt.xlabel('Real')
plt.ylabel('Complex')
plt.show()

Di seguito, l’output del programma: