Diary

The 1-D Linear Convection equation is the simplest, most basic model that can be used to learn something about CFD.

\(\frac{{\partial u}}{{\partial t}} + c\frac{{\partial u}}{{\partial x}} = 0\)

With given initial conditions (understood as a wave), the equation represents the propagation of that initial wave with speed \(c\), without change of shape.

Let the initial condition be \(u(x,0)=u_0(x)\). Then the exact solution of the equation is \(u(x,t)=u_0(x-ct)\).

We discretize this equation, in both space and time, using the Forward Difference scheme for the time derivative and the Backward Difference scheme for the space derivative.

Our discrete equation is:

\(\frac{{u_i^{n + 1} – u_i^n}}{{\Delta t}} + c\frac{{u_i^n – u_{i – 1}^n}}{{\Delta x}} = 0\)

Now let’s try implementing this in Python.

We’ll start by importing a few libraries to help us out.

• numpy is a library that provides a bunch of useful matrix operations akin to MATLAB;
• matplotlib is a 2D plotting library that we will use to plot our results;
• time and sys provide basic timing functions that we’ll use to slow down animations for viewing.

# Remember: comments in python are denoted by the pound sign
import numpy as np                      #here we load numpy
from matplotlib import pyplot as plt      #here we load matplotlib
import time, sys                   #and load some utilities



nx = 41         # spatial discretization: number of grid points
dx = 2 / (nx-1) # the distance between any pair of adjacent grid points
nt = 25         # nt is the number of timesteps we want to calculate
dt = .025       # dt is the amount of time each timestep covers (delta t)
c = 1           # assume wavespeed of c = 1

We have define an evenly spaced grid of points within a spatial domain that is 2 units of length wide, i.e., \(x_i \in [0,2]\). We have define a variable nx, which represent the number of grid points we want and dx is the distance between any pair of adjacent grid points.

We also need to set up our initial conditions. The initial velocity \(u_0\) is given as \(u=2\) in the interval \(0.5 \leq x \leq 1\) and \(u=1\) everywhere else in \((0,2)\) (i.e., a hat function). Here, we use the function ones() defining a numpy array which is nx elements long with every value equal to 1.

u = np.ones(nx)                      #numpy function ones()
u[int(.5 / dx):int(1 / dx + 1)] = 2  #setting u = 2 between 0.5 and 1 as per our I.C.s

x = np.linspace(0,2,nx)
plt.plot(x,u, 'red', label='Initial shape')

plt.show()

Now it’s time to implement the discretization of the convection equation using a finite-difference scheme.

We’ll store the result in a new (temporary) array un, which will be the solution \(u\) for the next time-step. We will repeat this operation for as many time-steps as we specify and then we can see how far the wave has convected.

We first initialize our placeholder array un to hold the values we calculate for the \(n+1\) timestep, using once again the NumPy function ones().

Then, we may think we have two iterative operations: one in space and one in time (we’ll learn differently later), so we’ll start by nesting one loop inside the other. Note the use of the nifty range() function. When we write: for i in range(1,nx) we will iterate through the u array, but we’ll be skipping the first element (the zero-th element).

un = np.ones(nx) #initialize a temporary array

for n in range(nt):        #loop for values of n from 0 to nt, so it will run nt times
    un = u.copy()          ##copy the existing values of u into un
    for i in range(1, nx): ## you can try commenting this line and...
        u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
    

plt.plot(x,u, 'green', label='Final shape')

plt.xlabel('space')
plt.ylabel('value')
pplt.title('1-D Linear Convection : nx=41, nt=25, dt=.025 and c=1')


plt.legend()

plt.show()

OK! So our hat function has definitely moved to the right, but it’s no longer a hat. What’s going on?

Note that:

\( – \left| x \right| \leq x \leq \left| x \right|\)

\( – \left| y \right| \leq y \leq \left| y \right|\)

After adding,

\( – \left( {\left| x \right| + \left| y \right|} \right) \leq x + y \leq \left( {\left| x \right| + \left| y \right|} \right)\)

Use the fact that \(\left| b \right| \leq a \Leftrightarrow – a \leq b \leq a \) (with \(b\) replaced by \(x+y\) and \(a\) by \(\left| x \right| + \left| y \right|\)), we have

\(\left| {x + y} \right| \leq \left| x \right| + \left| y \right|\)

Now, note that

\(\left| {x – y} \right| = \left| {\left( {x – z} \right) + \left( {z – y} \right)} \right|\)

Now we use tre triangle inequality and the fact that \(\left| {z – y} \right| = \left| {y – z} \right|\):

Estimate:

\(\int\limits_0^{\frac{\pi }{2}} {\frac{{{d} x}}{{1 + {{\tan }^{\sqrt 2 }}x}}}. \)

Solution: The problem cannot be evaluated by the usual techniques of integration; that is to say, the integrand does not have an antiderivative.

The problem can be handled if we happen to notice thet the integrand is symmetric about the point \(\left( {\frac{\pi }{4},\frac{1}{2}} \right)\).

To show this is so, let

\(f(x) = \frac{1}{{1 + {{\tan }^{\sqrt 2 }}x}}\)

It suffices to show that

\(f(x)+f(\pi /2 -x) = 1 \qquad \forall x \in \left[0, \pi/2\right]\)

Note: The diagram of \(f(\pi / 2 – x)\) is shown below:

It follow, from the symmetry just proved, that the area under the curve in \(\left[0, \pi/2\right]\) is one-half area of rectangle, i.e. \(\frac{1}{2} \frac{\pi}{2}= \frac{\pi}{4}\). So:

\(\int\limits_0^{\pi /2} {\frac{{dx}}{{1 + {{\tan }^{\sqrt 2 }}x}}} = \frac{\pi }{4}\)

In this post, we will deal with data from ECDC and we will explain basic data manipulation in Python with the Pandas package.

In our day, data is everywhere in enormous size and depth. Data science is an emerging field that penetrates every aspect of our life and, lately, it has proved to be an extraordinary weapon for predicting infections from Covid-19 and organizing strategies to limit the damage.

To import Pandas and Matplotlib packages we code:

import pandas as pd
import matplotlib.pyplot as plt

We download the excel file locally from ECDC site and open it using the read_excel function of Pandas library. We have named the file as data.xls in our case.

df=pd.read_excel("data.xlsx", engine="openpyxl")

We can first explore the data and the columns of the dataframe df:

We observe the columns of the dataframe — in our case, we will use the columns: dateRep, cases and deaths. Additionally, the name of the country is stored in column countriesAndTerritories.

We next select ‘Italy’ as the country under study. A new column is created named DateTime of type datetime where we store the day. In the following, we create a new dataframe with the name df_italia which is the same as the dataframe df_italia_sorted but it is sorted according to the column DateTime

df_italia=df[df.countriesAndTerritories=='Italy']
df_italia['DateTime']=pd.to_datetime(df_italia['dateRep'],format="%d/%m/%Y")

#We sort according to DateTime
df_italia_sorted=df_italia.sort_values(by='DateTime')

df_italia_selected=df_italia_sorted[df_italia_sorted.month>4]

We are interested in data after the month of April (i.e., May, June, July, August, … etc) so we choose to filter using the column month and create a new dataframe df_italia_selected.

Since the data in columns cases and deaths may have great variation, it is practical in order to understand the trend to use a moving average. We choose a moving average of seven days and we create two new columns (Moving Average Cases and Moving Average Deaths) where we store the average values of cases and deaths.

#Calculate moving average

df_italia_selected['Moving Average Cases']=df_italia_selected.cases.rolling(7,min_periods=1).mean()
df_italia_selected['Moving Average Deaths']=df_italia_selected.deaths.rolling(7,min_periods=1).mean()

We now plot the cases and deaths as functions of time. We choose the red color for cases and blue for deaths. It is useful to plot cases and deaths in the same figure with common x-axis in order to understand possible connection and relation. So, we use the subplots function and first create figure fig and axis ax1 (this will be the axis for the cases and it will be the left axis). We then create ax2 using twinx function. The values for deaths will be our right axis. A dashed line is used for the average values.

#Figure

fig, ax1=plt.subplots()
color1='tab:red'

ax1.plot(df_italia_selected['DateTime'], df_italia_selected['cases'], color=color1)

ax1.plot(df_italia_selected['DateTime'], df_italia_selected['Moving Average Cases'], color=color1,linestyle='dashed')

ax1.set_xlabel('Data')

ax1.set_ylabel('Cases',color=color1)
ax1.tick_params(axis='y',labelcolor=color1)

locs, labels=plt.xticks()
plt.setp(labels,rotation=90)

ax2=ax1.twinx() #instantiate a second axes that shares the same x-axis
color2='tab:blue'

ax2.plot(df_italia_selected['DateTime'], df_italia_selected['deaths'], color=color2)
ax2.plot(df_italia_selected['DateTime'], df_italia_selected['Moving Average Deaths'], color=color2,linestyle='dashed')

ax2.set_ylabel('Deaths',color=color2)
ax2.tick_params(axis='y',labelcolor=color2)

fig.tight_layout() #otherwise the right y-label is slightly clipped

The figure below is the program output.

^{Cases and deaths as a function of data for Italy}

Ci sono giorni in cui ogni cosa che vedo mi sembra carica di significati: messaggi che mi sarebbe difficile comunicare ad altri, definire, tradurre in parole, ma che appunto perciò mi si presentano come decisivi. Sono annunci o presagi che riguardano me e il mondo insieme: e di me non gli avvenimenti esteriori dell’esistenza ma ciò che accade dentro, nel fondo; e del mondo non qualche fatto particolare ma il modo d’essere generale di tutto. Comprenderete dunque la mia difficoltà a parlarne, se non per accenni.

Italo Calvino, da Se una notte d’inverno un viaggiatore

A prime number is a natural number greater than 1 that is not a product of two smaller natural numbers.

Euclid proved the infinity of primes by contradiction.

• Assume there are a finite number, n , of primes , the largest being \(p_n\);

• Consider the number that is the product of these, plus one: \(N = \prod\limits_{i = 1}^n {{p_i}} + 1\)

• By construction, \(N\) is not divisible by any of the \(p_i\).

• Hence it is either prime itself, or divisible by another prime greater than \(p_n\), contradicting the assumption.

Euclid’s proof is the easiest to understand, especially if you’re not well-versed in mathematics, but now I am going to introduce another proof of the infinitude of primes, which is shorter and I also find it to be the most appealing, at least in my eyes.

In mathematics, the prime-counting function is the function counting the number of prime numbers less than or equal to some real number \(x\). E.g. for \(x=10.124\), we have that 4 number of prime number: 2, 3, 5, 7.

Let \(\pi(x)\) be a prime-counting function.

The Prime Number Theorem suggests that:

\(\mathop {\lim }\limits_{x \to \infty } \frac{{\pi \left( x \right)}}{{x/\ln (x)}} = 1\)

In other words, we can say that \(\pi(x)\) “behaves” like \(x/\ln (x)\) when \(x\) is approaching infinity — there’s an asymptotic ‘equality’.

So, there are infinite primes.

---

For Python code, we use SymPy – a Python library for symbolic mathematics. In particular, we’ll use primerange(a, b) for generate all prime numbers in the range [a, b).

The code:

import sympy
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as ani
#plt.style.use('dark_background')
plt.style.use('default')

def convert_to_polar_coordinates(num):
    return num*np.cos(num), num*np.sin(num)


lim = 10000

fig, ax = plt.subplots(figsize = (8, 8))
primes_numbers = sympy.primerange(0, lim)
primes_numbers = np.array(list(primes_numbers))
x, y = convert_to_polar_coordinates(primes_numbers)


def init():
    ax.cla()
    ax.plot([], [])


def plot_in_polar(i):
    ax.cla()
    ax.plot(x[:i], y[:i], linestyle='', marker='o', markersize=0.75, c='#FC0000')
    ax.set_xlim(-lim, lim)
    ax.set_ylim(-lim, lim)
    ax.axis("off")

plt.savefig('prime.png')

animator = ani.FuncAnimation(fig=fig, func=plot_in_polar, interval=1, frames=len(primes_numbers))

animator.save("prime.gif")

plt.show()

Sia \(d\) un numero e \(n \in {\mathbb{N}^ + }\) tale che:

\({n^2} < d < {\left( {n + 1} \right)^2}\)

ovvero:

\(n < \sqrt d < \left( {n + 1} \right)\)

e quindi:

\(0< \sqrt{d}-n<1\)

Per assurdo, sia

\(\sqrt d = \frac{p}{q}{\text{ con }}p,q \in {\mathbb{N}^ + }\)

Più precisamente, assumeremo che \(q\) sia il minimo dei possibili denominatori per queste espressioni razionali di \(\sqrt d\) e quindi assumiamo che sia il minimo \(q\) tale che \(q \sqrt{d}=p \in \mathbb{N}^{+}\).

Quindi:

\(\left( {\sqrt d – n} \right)q\sqrt d = qd – nq\sqrt d \in {\mathbb{N}^ + }\)

e pertanto, essendo \(\sqrt{d}-n<1\) dev’essere:

\(\left( {\sqrt d – n} \right)q < q\)

contro il fatto che, per ipotesi, \(q\) è il minimo dei possibili denominatori per la rappresentazione razionale di \(\sqrt d\).