text:
```python
import theme
theme.load_style()
```

# Lesson 4: Plane and Space Trusses

This lecture by Tim Fuller is licensed under the
Creative Commons Attribution 4.0 International License.  All code examples are also licensed under the [MIT license](http://opensource.org/licenses/MIT).

# Topics

- [Introduction](#intro)
- [Local and Global Coordinates](#loc_glob_coords)
  - [Overview](#lg_overview)
  - [Coordinate Transformation](#coord_trans)
  - [Truss Equations in Global Coordinates](#g_trusseq)
- [Example](#ex1)
  - [Element Connectivity](#el_conn_ex1)
  - [Boundary Conditions](#bcs_ex1)
  - [Solution](#sol_ex1)
- [Computational Implementation](#computational_implementation)

# Introduction[](#top)

Are the above structures frames or trusses?

**Frame**

- Carries transverse and bending loads
- Members are connected rigidly

**Truss**

- Members only carry axial loads (no transverse shearing or bending)
- Members are connected by pins that allow free rotation

Because of their pin connections, individual members of truesses support only axial loads, as shown below.

Because of this property, members of trusses are termed "two force" members.  Importantly, though members in a truss support only loads along their axis, the truss is capable of supporting loads in the $XY$ plane (2D plane trusses) or $XYZ$ space (3D space trusses).  This gives rise to the notion of a "local" coordinate system for individual members and a a "global" coordinate system for the truss.  In general, force-displacement relationships for members are formulated in their local coordinate system and transformed to the global truss coordinate system where they are solved.  The remainder of this chapter will focus on formulating and solving these systems of equations for trusses.

# Local and Global Coordinate Systems[](#top)

## Overview

A member of a truss is a 1D two-force member embedded in a 2D (plane truss) or 3D (space truss) coordinate system, as shown below.  In their local coordinates ($\hat{x}$, $\hat{y}$), members of a truss experience forces and displacements in only their local $\hat{x}$ direction.  When viewed in the global coordinates, members experience forces and displacements in the $x$ and the $y$ coordinate directions.

To accomodate the extra dimensionality, rows and columns of zeros are inserted in the local element stiffness equation corresponding to the (non-existent) forces and displacements in the local $\hat{y}$ direction.

1D element coordinate system

  $$
  \begin{Bmatrix} f_{i\hat{x}} \\ f_{j\hat{x}} \end{Bmatrix} =
  k^{(e)}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
  \begin{Bmatrix} u_{i\hat{x}} \\ u_{j\hat{x}} \end{Bmatrix}
  $$

2D element coordinate system

  $$
  \begin{Bmatrix}
    f_{i\hat{x}} \\ f_{i\hat{y}} \\ f_{j\hat{x}} \\ f_{j\hat{y}}
  \end{Bmatrix} =
  k^{(e)}\begin{bmatrix}
    1 & 0 & -1 & 0 \\
    0 & 0 & 0 & 0 \\
    -1 & 0 & 1 & 0 \\
     0 & 0 & 0 & 0 \\
  \end{bmatrix}
  \begin{Bmatrix}
    u_{i\hat{x}} \\ u_{i\hat{y}} \\ u_{j\hat{x}} \\ u_{j\hat{y}}
  \end{Bmatrix}$$

Casting the local equations in the global 2D system requires a change of coordinates by means of [Coordinate Transformations](#coord_trans).

## Coordinate Transformation

Consider the vector $\boldsymbol{v}$ represented in a local coordinate system ($\hat{x},\hat{y}$) embedded in a global coordinate system ($x,y$).

The components of $\boldsymbol{v}$ in the local coordinate system $v_{\hat{x}}$ and $v_{\hat{y}}$ are related to the components $v_x$ and $v_y$ in the global coordinates by

$$
\begin{align}
v_x &= v_{\hat{x}} \cos \theta - v_{\hat{y}} \sin \theta \\
v_y &= v_{\hat{x}} \sin \theta + v_{\hat{y}} \cos \theta
\end{align}
$$

In matrix form, we say

$$
\begin{Bmatrix} v_x \\ v_y \end{Bmatrix} =
\begin{bmatrix}
\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta
\end{bmatrix}
\begin{Bmatrix} v_{\hat{x}} \\ v_{\hat{y}} \end{Bmatrix}
$$

we recognize

$$
R_{ij} = \begin{bmatrix}
\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta
\end{bmatrix}
$$

as the familiar rotation matrix.

The components of displacement and force in the local coordinates transform similarly.  In compact notation we right

Displacement relationship

$$
\begin{align}
u_{ix} &= u_{i\hat{x}} \cos \theta - u_{i\hat{y}} \sin \theta \\
u_{iy} &= u_{i\hat{x}} \sin \theta + u_{i\hat{y}} \cos \theta \\
u_{jx} &= u_{j\hat{x}} \cos \theta - u_{j\hat{y}} \sin \theta \\
u_{jy} &= u_{j\hat{x}} \sin \theta + u_{j\hat{y}} \cos \theta
\end{align}
$$

Force relationship

$$
\begin{align}
f_{ix} &= f_{i\hat{x}} \cos \theta - f_{i\hat{y}} \sin \theta \\
f_{iy} &= f_{i\hat{x}} \sin \theta + f_{i\hat{y}} \cos \theta \\
f_{jx} &= f_{j\hat{x}} \cos \theta - f_{j\hat{y}} \sin \theta \\
f_{jy} &= f_{j\hat{x}} \sin \theta + f_{j\hat{y}} \cos \theta
\end{align}
$$

In matrix notation, we write

$$
\begin{Bmatrix} u_{ix} \\ u_{iy} \\ u_{jx} \\ u_{jy} \end{Bmatrix} =
\begin{bmatrix}
\cos \theta & -\sin \theta & 0 & 0 \\
\sin \theta & \cos \theta & 0 & 0 \\
0 & 0 & \cos \theta & -\sin \theta \\
0 & 0 & \sin \theta & \cos \theta
\end{bmatrix}
\begin{Bmatrix} u_{i\hat{x}} \\ u_{i\hat{y}} \\ u_{j\hat{x}} \\ u_{j\hat{y}} \end{Bmatrix}
$$

and similarly for the force.

Remark

Each node now has 2 degrees of freedom

 Remark

Note the ordering of the nodal displacements for the now multi-DOF displacement array above, repeated here

$$
\begin{Bmatrix}
u_{1x} \\
u_{1y} \\
u_{2x} \\
u_{2y} \\
\vdots \\
u_{nx} \\
u_{ny} \\
\end{Bmatrix}
\begin{matrix}
\longleftarrow 1^{\text{st}} \text{ and } 2^{\text{nd}} \text{ DOF for node 1} \\[8pt]
\longleftarrow 1^{\text{st}} \text{ and } 2^{\text{nd}} \text{ DOF for node 2} \\[7pt]
\vdots \\[7pt]
\longleftarrow 1^{\text{st}} \text{ and } 2^{\text{nd}} \text{ DOF for node n}
\end{matrix}
$$

Other node ordering schemes exit, such as

$$
\begin{Bmatrix}
u_{1x} \\
u_{2x} \\
\vdots \\
u_{nx} \\
u_{1y} \\
u_{2y} \\
\vdots \\
u_{ny} \\
\end{Bmatrix}
$$

Each ordering technique has its merit.  In this course, we adopt the former ordering method.

With the chosen array ordering, the global degree of freedom $DOF$ of the local degree of freedom $dof$ of node $i$ is given by

$$DOF = 2i + dof$$

## Truss Equations in the Global Coordinates[](#top)

Returning to the force displacement relationships.  The transformation from local to global coordinates for the displacement and force arrays is given by

$$
\{u\} = [T]\{\hat{u}\}
$$
$$
\{f\} = [T]\{\hat{f}\}
$$

where

$$
T = \begin{bmatrix}
\cos \theta & -\sin \theta & 0 & 0 \\
\sin \theta & \cos \theta & 0& 0 \\
0 & 0 & \cos \theta & -\sin \theta \\
0 & 0 & \sin \theta & \cos \theta
\end{bmatrix}
$$

In the local coordinates, recall that

$$
\{\hat{f}\} = [\hat{k}] \{\hat{u}\}
$$

or,

$$
\{\hat{f}\} = [T]^{-1} \{f\} = [\hat{k}] \{\hat{u}\} = [\hat{k}] [T]^{-1} \{u
\}$$

which simplifies to

$$
\{f\} = [T] \{\hat{u}\} [T]^{-1} \{u\}
$$

or

$$
\{f\} = [k] \{u\}
$$

where

$$
[k] = [T] \{\hat{k}\} [T]^{-1}
$$

You'll recall from linear algebra that $[T]$ is orthogonal, so $[T]^{-1} = [T]^{T}$

Following through on the multiplication leads to

$$
[k] = \hat{k}\begin{bmatrix} \cos^{2} \theta & \sin \theta \cos \theta  & -\cos^{2} \theta & -\sin \theta \cos \theta \\ \sin \theta \cos \theta & \sin^{2} \theta & - \sin \theta \cos \theta & - \sin^{2} \theta \\ - \cos^{2} \theta & - \sin \theta \cos \theta & \cos^{2} \theta & \sin \theta \cos \theta \\ - \sin \theta \cos \theta & - \sin^{2} \theta & \sin \theta \cos \theta & \sin^{2} \theta\end{bmatrix}
$$

It turns out that $[k]$ has the same structure as $[\hat{k}]$, namely

$$
[k] = \hat{k}\begin{bmatrix} a & b \\ -b & a \end{bmatrix}
$$

In $[\hat{k}]$, $a=b=1$ whereas in $[k]$, $a$ and $b$ are matrices.  Being that the two views of the element stiffness have the same structure, the assembly of the components in to the global stiffness also has the same structure, as shown in the images below for five element structures in 1 and 2D.

Note that the assembly procedure is the same and that, generally speaking, the shape of $[K]$ is $n_{\text{dof}}^{\text{sys}} \times n_{\text{dof}}^{\text{sys}}$, $n_{\text{dof}}^{\text{sys}} = n\cdot n_{\text{dof}}^{\text{n}}$ where $n$ and $n_{\text{dof}}^{\text{n}}$ are the number of nodes and the number of degrees of freedom per node, respectively.  We see then that the assembly procedure in 1D is just a special case of the assembly procedure in 2D (any D, really).  We'll see how to exploit this fact to generalize our existing finite element solver to solve systems in multiple dimensions.

### Alternative Method of Assembling $[k]$

It turns out that we need to explicitly compute all the $\theta$ to compute $[k]$.  Recall

So,

$$
[k] = [T] [\hat{k}] [T]^{T} = \hat{k}\begin{bmatrix}
n_{x}^{2} & n_{x} n_{y} & -n_{x}^{2} & -n_{x} n_{y}\\
n_{x} n_{y} & n_{y}^{2} & -n_{x} n_{y} & -n_{y}^{2}\\
-n_{x}^{2} & -n_{x} n_{y} & n_{x}^{2} & n_{x} n_{y}\\
-n_{x} n_{y} & -n_{y}^{2} & n_{x} n_{y} & n_{y}^{2}
\end{bmatrix}
$$

The matrix above can be computed from the *outer product* of $\boldsymbol{n}$ with itself.  The outer product (often called the dyadic product) is denoted in some texts as $\boldsymbol{a} \otimes \boldsymbol{b}$ and is defined as

$$
\boldsymbol{C} = \boldsymbol{a} \otimes \boldsymbol{b}
\qquad C_{ij} = a_i b_j
$$

#### Calculating Element Unit Vectors

Exactly how do we calculate the unit vectors needed by the previous relationship?  Consider the following truss, in particular element 3

The nodes $i$ and $j$ associated with each element come from the element connectivity table.  $x_i$ and $x_j$ come from the table of nodal coordinates.

The element direction vector is given by

$$
v^{(3)} = x_j - x_i
$$

and the element length is

$$
l^{(3)} = || v^{(3)} || = \sqrt{v^{(3)} \cdot v^{(3)}} = \sqrt{v_x^2 + v_y^2}
$$

Finally, the element unit vector is

$$
n^{(3)} = \frac{v^{(3)}}{l^{(3)}}
$$

# Example[](#top)

Let's put it together with an example.  Let us consider the following truss.  We begin as with the 1D problems by analyzing the element connectivity and determining boundary conditions

## Element Connectivity

Assembling the stiffness matrix is analogous to assembling 1D systems.  We use a connectivity table as before

| Element Number  | Node i | Node j |
| :-------------: |:------:| :-----:|
| 1               | 1      | 2      |
| 2               | 1      | 3      |

##  Boundary Conditions

### Pinned Boundary Conditions

The boundary conditions on the nodes 2 and 3 are examples of pinned boundary conditions - not allowing moving in the $x$ or the $y$ directions.  Or, $u_{2x} = u_{2y} = u_{3x} = u_{3y} = 0$.

### Roller Boundary Conditions

The boundary condition on $x$ dof of node 1 is an example of a roller boundary condition.  The displacement $u_{1y}$ is a model *unknown*

### Displacement Boundary Conditions

The boundary condition on the $y$ dof of noe 1 is an example of a displacement boundary condition.  The displacement $u_{1x} = -\delta$.

### Force Boundary Conditions

A nodal force is applied at node 1: $f_{1y} = P$

##  Solution

The solution to the above problem is worked out in the computational cell below.  One new thing to note is the use of the `femlib.viewer.display_truss` method that produces plots of truss members.  See the usage below and the `femlib.viewer.display_truss` docstring for details.  Also note the use of `-1` in the call to `np.reshape`.  See the `numpy` [reshape](http://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html) documentation for details.

```python
import numpy as np
from sympy import symbols, cos, sin, init_printing, Matrix
init_printing()
```

```python
# allows plots to be displayed inline
%matplotlib inline
from femlib.ifem import display_truss

num_node = 3
num_dof_per_node = 2
num_elem = 2
num_dof = num_node * num_dof_per_node

# vertices[i,j] is the jth coordinate of the ith node
node_map = [1, 2, 3]
vertices = np.array([[0., 0.],
                     [3., 4.],
                     [0., 4.]])

# connect[i,j] is the jth node number of the ith element
connect = np.array([[1, 2], [1, 3]])

# material properties
E, A = 70E+09, 5 * .01 * .01

# Load
P = 1000E+03

# Boundary conditions
delta = .05

# Element 1 stiffness
v1 = vertices[1] - vertices[0]
l1 = np.sqrt(np.dot(v1, v1))
n1 = v1 / l1
n1n1 = np.outer(n1, n1)
k1 = A * E / l1 * n1n1

# Element 2 stiffness
v2 = vertices[2] - vertices[0]
l2 = np.sqrt(np.dot(v2, v2))
n2 = v2 / l2
n2n2 = np.outer(n2, n2)
k2 = A * E / l2 * n2n2

# Global stiffness
K = np.zeros((num_dof, num_dof))

# Element 1 contribution to global stiffness.  Note, that
# K1 is a 2x2 so that we are adding in the full 4x4 element
# stiffness to the global stiffness
i, j = connect[0] - 1
I, J = num_dof_per_node * i, num_dof_per_node * j
K[I:I+2, I:I+2] += k1
K[J:J+2, J:J+2] += k1
K[I:I+2, J:J+2] -= k1
K[J:J+2, I:I+2] -= k1

# Element 2 contribution to global stiffness
i, j = connect[1] - 1
I, J = num_dof_per_node * i, num_dof_per_node * j
K[I:I+2, I:I+2] += k2
K[J:J+2, J:J+2] += k2
K[I:I+2, J:J+2] -= k2
K[J:J+2, I:I+2] -= k2

# Global force
F = np.zeros(num_dof)

# P applied to node 1 in y direction
node, dof = 0, 1
F[num_dof_per_node * node + dof] = P

# Apply boundary conditions via penalty method
X = 1.e30
Kbc = K.copy()
Fbc = F.copy()

# node 1: -delta applied in x direction
Kbc[0,0] = X
Fbc[0] = -delta * X

# nodes 2,3: 0 applied to x,y directions
Kbc[2,2] = Kbc[3,3] = Kbc[4,4] = Kbc[5,5] = X
Fbc[2] = Fbc[3] = Fbc[4] = Fbc[5] = 0 *X

# Solve for displacement
u = np.linalg.solve(Kbc, Fbc)

# Find reactions
r = np.dot(K, u) - F

# Updated nodal positions
x = vertices + u.reshape((num_node, -1))

# Display the truss
display_truss(vertices, connect)

# display the results
display_truss(vertices, connect, overlay=x, oc="b", ol="Deformed")

print u.reshape((num_node, -1))
```

# Computational Implementation[](#top)

Would it be possible to modify `oned.solve_system` in such a way that it would be able to solve a truss problem and still be compatible with our `uniform_bar` and `tapered_bar` applications?  Yes!  And it requires only a slight modification to the assembly procedures to boot.

Let's investigate the current `oned.solve_system` code and identify places where it will need to be modified, starting with the assembly of the global stiffness

K = np.zeros((num_dof, num_dof))
for (el, nodes) in enumerate(connect):
    i = dof_map(nodes[0])
    j = dof_map(nodes[1])
    K[i, i] += k[el]
    K[j, j] += k[el]
    K[i, j] -= k[el]
    K[j, i] -= k[el]

 Note

The implicit assumption that k[el] is a scalar is made in the algorithm.  For problems in more than 1 dimension (like a plane truss), k[el] is instead a matrix.

The question then arises: how do we modify the above algorithm to allow for `k[el]` to be a matrix?  The answer is quite simple, change the left hand side of the assignment, e.g, `K[i, i] += k[el]` to use the appropriate slice of `K` and not just a single index.  The details are left as an exercise at the end of the chapter.

It turns out the modifying the global stiffness assembly is the only change needed to convert our existing code to work with trusses!

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: pinkieli/fem-with-python

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: pinkieli/fem-with-python

max_forks_repo_licenses: MIT Unlicense


text:
> Texto fornecido sob a Creative Commons Attribution license, CC-BY. Todo o código está disponível sob a FSF-approved BSD-3 license.

> (c) Original por Lorena A. Barba, Gilbert F. Forsyth em 2017, traduzido por Felipe N. Schuch em 2020.

> [@LorenaABarba](https://twitter.com/LorenaABarba) - [@fschuch](https://twitter.com/fschuch)

12 passos para Navier-Stokes
======
***

Continuamos nossa jornada para resolver a equação de Navier-Stokes, esse é o Passo 4. Não continue a menos que você tenha completado as tarefas anteriores. Na verdade, esse passo será uma combinação dos dois anteriores \[[2](./02_Passo_2.ipynb),[3](./04_Passo_3.ipynb)\]. O maravilhoso *reaproveitamento de código*!

Passo 4: Equação de Burgers
----
***

Você pode ler mais sobre a equação de Burgers na sua [página wikipedia](http://en.wikipedia.org/wiki/Burgers'_equation).

A equação de Burgers em uma dimensão se parece com:

$$\frac{\partial u(x,t)}{\partial t} + u(x,t) \frac{\partial u(x,t)}{\partial x} = \nu \frac{\partial ^2u(x,t)}{\partial x^2}$$

Como você pode perceber, ela é uma combinação entre convecção não linear e difusão. É surpreendente o quanto você aprende com essa pequena e elegante equação!

Discretizamos a equação usando os mesmos métodos que já foram detalhados nos Passos [1](./01_Passo_1.ipynb) a [3](./04_Passo_3.ipynb). Usando diferença para frente para a derivada temporal, diferença para trás para a derivada espacial de primeira ordem e um esquema de diferença centrada para a derivada segunda. A equação é então escrita como:

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} + u_i^n \frac{u_i^n - u_{i-1}^n}{\Delta x} = \nu \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}$$

Como antes, uma vez conhecida a condição inicial, a única incógnita é $u_i^{n+1}$. Vamos prosseguir no tempo com:

$$u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$$

### Condição Inicial e de Contorno

Para examinar algumas propriedades interessantes da equação de Burgers, é adequado o uso de condições inicial e de contorno diferentes do que usamos nas tarefas anteriores.

A nossa condição inicial para esse problema agora será:

\begin{eqnarray}
u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\
\phi &=& \exp \bigg(\frac{-x^2}{4 \nu} \bigg) + \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu} \bigg)
\end{eqnarray}

Que tem solução analítica, dada por:

\begin{eqnarray}
u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\
\phi &=& \exp \bigg(\frac{-(x-4t)^2}{4 \nu (t+1)} \bigg) + \exp \bigg(\frac{-(x-4t -2 \pi)^2}{4 \nu(t+1)} \bigg)
\end{eqnarray}

Nossa condição de contorno será:

$$u(0) = u(2\pi)$$

Ela é denominada condição de contorno *periódica*. Preste atenção! Isso pode lhe causar certa dor de cabeça se você não tratá-la cuidadosamente.

### Poupando tempo com SymPy

A condição inicial que usaremos para a Equação de Burgers pode ser meio trabalhosa para ser calculada à mão. A derivada $\frac{\partial \phi}{\partial x}$ não é lá muito complicada, mas seria fácil perder um sinal ou esquecer um fator de $x$ em algum lugar, então vamos usar Sympy para nos ajudar.

[SymPy](http://sympy.org/en/) é uma biblioteca Python para metemática simbólica. Ela tem várias das funcionalidades de cálculo simbólico como Mathematica, com o benefício adicional de que podemos facilmente converter os resultados de volta para nosso cálculo em Python (SymPy também é grátis e de código aberto).

Começamos carregando a biblioteca SymPy, em conjunto com nossa biblioteca preferida, Numpy.

```python
import numpy
import sympy
```

Vamos também informar ao SymPy que nós queremos que todas impressões sejam renderizadas usando $\LaTeX$. Isso vai deixar nosso Notebook incrível!

```python
from sympy import init_printing
init_printing(use_latex=True)
```

Começamos por definir nossas variáveis simbólicas para as três variáveis presentes na nossa condição inicial e então digitar a equação completa para $\phi$. Vamos obter uma versão satisfatoriamente renderizada para nossa equação para $\phi$.

```python
x, nu, t = sympy.symbols('x nu t')
phi = (sympy.exp(-(x - 4 * t)**2 / (4 * nu * (t + 1))) +
       sympy.exp(-(x - 4 * t - 2 * sympy.pi)**2 / (4 * nu * (t + 1))))
phi
```

Ela parece um pouco pequena, mas está tudo em ordem. Agora, obter a derivada parcial $\frac{\partial \phi}{\partial x}$ é uma tarefa trivial.

```python
phiprime = phi.diff(x)
phiprime
```

Se você quiser ver a versão não rederizada, basta apenas usar o comando `print()`.

```python
print(phiprime)
```

    -(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))

### E agora?

Uma vez que temos a versão *Pythonica* da nossa derivada, podemos terminar de escrever a condição inicial completa e então traduzir ela para uma expressão Python utilizável. Para tanto, usaremos a função *lambdify*, que toma uma equação simbólica do SymPy e a transforma em uma função invocável.

```python
from sympy.utilities.lambdify import lambdify

u = -2 * nu * (phiprime / phi) + 4
print(u)
```

    -2*nu*(-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)))/(exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1))) + exp(-(-4*t + x)**2/(4*nu*(t + 1)))) + 4

### Lambdify

Para que lambdify transforme a expressão em uma função invocável, dizemos ao lambdify quais variáveis deve solicitar e a equação na qual queremos conectá-la.

```python
ufunc = lambdify((t, x, nu), u)
print(ufunc(1, 4, 3))
```

    3.49170664206445

### De volta à equação de Burgers

Agora que temos a condição inicial preparada, podemos prosseguir para finalizar nosso problema. Podemos exibir nossa condição inicial no gráfico da condição inicial (CI) usando nossa função lambdify.

```python
from matplotlib import pyplot
%matplotlib inline

###Declaração de variáveis
nx = 101
nt = 100
dx = 2. * numpy.pi / (nx - 1)
nu = .07
dt = dx * nu

x = numpy.linspace(0, 2 * numpy.pi, nx)
un = numpy.empty(nx)
t = 0

u = numpy.asarray([ufunc(t, x0, nu) for x0 in x])
u
```

    array([4.        , 4.06283185, 4.12566371, 4.18849556, 4.25132741,
           4.31415927, 4.37699112, 4.43982297, 4.50265482, 4.56548668,
           4.62831853, 4.69115038, 4.75398224, 4.81681409, 4.87964594,
           4.9424778 , 5.00530965, 5.0681415 , 5.13097336, 5.19380521,
           5.25663706, 5.31946891, 5.38230077, 5.44513262, 5.50796447,
           5.57079633, 5.63362818, 5.69646003, 5.75929189, 5.82212374,
           5.88495559, 5.94778745, 6.0106193 , 6.07345115, 6.136283  ,
           6.19911486, 6.26194671, 6.32477856, 6.38761042, 6.45044227,
           6.51327412, 6.57610598, 6.63893783, 6.70176967, 6.76460125,
           6.82742866, 6.89018589, 6.95176632, 6.99367964, 6.72527549,
           4.        , 1.27472451, 1.00632036, 1.04823368, 1.10981411,
           1.17257134, 1.23539875, 1.29823033, 1.36106217, 1.42389402,
           1.48672588, 1.54955773, 1.61238958, 1.67522144, 1.73805329,
           1.80088514, 1.863717  , 1.92654885, 1.9893807 , 2.05221255,
           2.11504441, 2.17787626, 2.24070811, 2.30353997, 2.36637182,
           2.42920367, 2.49203553, 2.55486738, 2.61769923, 2.68053109,
           2.74336294, 2.80619479, 2.86902664, 2.9318585 , 2.99469035,
           3.0575222 , 3.12035406, 3.18318591, 3.24601776, 3.30884962,
           3.37168147, 3.43451332, 3.49734518, 3.56017703, 3.62300888,
           3.68584073, 3.74867259, 3.81150444, 3.87433629, 3.93716815,
           4.        ])

```python
pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.plot(x, u, marker='o', lw=2)
pyplot.xlim([0, 2 * numpy.pi])
pyplot.ylim([0, 10]);
```

Essa definitivamente não é a função chapéu com a qual viemos lidando até agora. Chamamo ela de *função dente de serra* ("saw-tooth function"). Seguimos adiante para ver o que acontece.

### Condição de Contorno Periódica

Uma das grandes diferenças entre o Passo 4 e as lições anteriores é o emprego da condição de contorno *periódica*. Se você experimentar diferentes parâmetros nos Passos 1 e 2 e tentar executar as simulações por mais tempo (aumentando `nt`) vai notar que a onda inicial vai continuar se movendo para a direita até que ela não mais sequer aparece no gráfico.

Com condições de contorno periódicas, quando um ponto chega ao lado direito do domínio, ele *contorna* de volta para o início do domínio.

Apressentamos novamente a equação discretizada que desenvolvemos no início desse notebook:

$$
u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)
$$

Qual o *significado* de $u_{i+1}^n$ quando $i$ já está ao final do domínio?

Pense sobre isso por um minuto antes de proceguir.

```python
for n in range(nt): #Laço temporal
    un = u.copy()
    for i in range(1, nx-1):
        u[i] = un[i] - un[i] * dt / dx *(un[i] - un[i-1]) + nu * dt / dx**2 *\
                (un[i+1] - 2 * un[i] + un[i-1])
    u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu * dt / dx**2 *\
                (un[1] - 2 * un[0] + un[-2])
    u[-1] = u[0]

u_analytical = numpy.asarray([ufunc(nt * dt, xi, nu) for xi in x])
```

```python
pyplot.figure(figsize=(11, 7), dpi=100)
pyplot.plot(x,u, marker='o', lw=2, label='Computacional')
pyplot.plot(x, u_analytical, label='Analítico')
pyplot.xlim([0, 2 * numpy.pi])
pyplot.ylim([0, 10])
pyplot.legend();
```

***

O que vem agora?
----

Os Passos que seguem, do 5 ao 12, serão considerados em duas dimensões. Mas será fácil extender os esquemas 1D de diferenças finitas para 2D ou 3D. Apenas aplique a definição: uma derivada parcial com respeiro a $x$ é a variação na direção $x$ *enquanto mantemos $y$ constante*.

Antes de proceguir para o [Passo 5](./07_Passo_5.ipynb), tenha certeza que tenha completado seu próprio código para os Passos 1 à 4, e que tenha experimentado diferentes combinações de parâmetros e refletido sobre o que está acontecendo. Adicionalmente, recomendamos que faça uma pequena pausa para aprender sobre [operações com arranjos com NumPy](./06_Operacoes_de_arranjos_com_NumPy.ipynb).

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>
<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>

> A célula acima executa o estilo para esse notebook. Nós modificamos o estilo encontrado no GitHub de [CamDavidsonPilon](https://github.com/CamDavidsonPilon), [@Cmrn_DP](https://twitter.com/cmrn_dp).

ext: ipynb

max_stars_repo_name: fschuch/CFDPython

max_issues_repo_name: fschuch/CFDPython

max_forks_repo_name: fschuch/CFDPython


text:
```python
import numpy as np
from sympy import *
from sympy.calculus.util import continuous_domain
import pandas as pd
import copy
from P1 import *

init_printing(use_unicode=True)

# a.B

# matrix  = Matrix([
#     [8.8,-1.3,-6.8, -36.4],
#     [9.7, 4.5, 5.6, 108],
#     [7.5,-0.5,-0.8,9.06],
# ])

matrix = Matrix([
[7.40, 6.00, 0.700, -42.1],
[-6.70, -5.20, 5.60, 22.2],
[-2.70, -1.40, 7.40, -6.47],
])

# matrix = Matrix([
# [5.10, 3.20, -7.3, -136 ],
# [-8.40, -5.20, 8.10, 186],
# [4, 7.10, 0.4, -74]
# ])
# matrix = Matrix([[0.600, -9.90, -2.20, -47.2], [-3.70, -0.500, 1.80, -1.96], [4.4,8.10,1.8, 20.5]])

# pprint(escalona_gauss(matrix, 3, True)[0])
a = gauss_elimi(matrix, 3, True, 0.001)
#a.recursi()

# vec = Matrix([-9.390888310, -5.671787784, 9.583116262])
# vec = Matrix([-3.51584937, 1.021725838, 0.4840877769])
vec = Matrix([-1.788673724, -4.532451158, -2.384439280])
# # for x in a.xs:
# #     pprint(x-vec)
abs(a.xs[0] - vec)
```

    Interacao  0
    Antes da iteracao
    ⎡7.4   6.0   0.7  -42.1⎤
    ⎢                      ⎥
    ⎢-6.7  -5.2  5.6  22.2 ⎥
    ⎢                      ⎥
    ⎣-2.7  -1.4  7.4  -6.47⎦
    Depois da iteracao
    matriz
    ⎡7.4  6.0   0.7   -42.1⎤
    ⎢                      ⎥
    ⎢ 0   0.23  6.23  -15.9⎥
    ⎢                      ⎥
    ⎣ 0   0.79  7.66  -21.9⎦
    []
    multiplicadores
    ⎡  0     0  0  0⎤
    ⎢               ⎥
    ⎢-0.905  0  0  0⎥
    ⎢               ⎥
    ⎣-0.365  0  0  0⎦

    Interacao  1
    Antes da iteracao
    ⎡7.4  6.0   0.7   -42.1⎤
    ⎢                      ⎥
    ⎢ 0   0.23  6.23  -15.9⎥
    ⎢                      ⎥
    ⎣ 0   0.79  7.66  -21.9⎦
    Depois da iteracao
    matriz
    ⎡7.4  6.0   0.7   -42.1⎤
    ⎢                      ⎥
    ⎢ 0   0.79  7.66  -21.9⎥
    ⎢                      ⎥
    ⎣ 0    0    4.0   -9.53⎦
    [(1, 2, 'Col 1')]
    multiplicadores
    ⎡  0       0    0  0⎤
    ⎢                   ⎥
    ⎢-0.365    0    0  0⎥
    ⎢                   ⎥
    ⎣-0.905  0.291  0  0⎦

    Interacao  2
    Antes da iteracao
    ⎡7.4  6.0   0.7   -42.1⎤
    ⎢                      ⎥
    ⎢ 0   0.79  7.66  -21.9⎥
    ⎢                      ⎥
    ⎣ 0    0    4.0   -9.53⎦
    Depois da iteracao
    matriz
    ⎡7.4  6.0   0.7   -42.1⎤
    ⎢                      ⎥
    ⎢ 0   0.79  7.66  -21.9⎥
    ⎢                      ⎥
    ⎣ 0    0    4.0   -9.53⎦
    [(1, 2, 'Col 1')]
    multiplicadores
    ⎡  0       0    0  0⎤
    ⎢                   ⎥
    ⎢-0.365    0    0  0⎥
    ⎢                   ⎥
    ⎣-0.905  0.291  0  0⎦

$\displaystyle \left[\begin{matrix}0.128673724\\0.147548842\\0.00443928000000016\end{matrix}\right]$

```python

max(abs(a.xs[0] - vec))
```

```python
a.xs[0]

```

$\displaystyle \left[\begin{matrix}-1.66\\-4.68\\-2.38\end{matrix}\right]$

```python
A = Matrix([
    [-0.3, -0.3, 6.5, 21.3],
    [3.6, 3.4, -1.6, -10.7],
    [7.8, -1, 1.5, -16.6]
])

Real = Matrix([-2.601908433, 1.117749365, 3.208423428])

a = gauss_elimi(A, 3, True, 0.1)
a.recursi()

print(max(abs(a.xs[0] - Real)))
```

    achou em 1 iteracoes
    0.00225063500000000

```python
# A = Matrix([
#     [-4, 9, -9, -6],
#     [1, 3, 8, -6],
#     [10, -8, 3, -5]
# ])

# A = Matrix([
# [6, -11, 5, -3],
# [-2, 7, -10, 10],
# [12, -2, -10, 7]
# ])
#
#
#
#
# erro = 0.00100
# var = 13.8
# Ms = []
#
# Aswa = (0,1)
# Bswa = (0,2)
# Cswa = (1,2)
#
# swaps = [Aswa, Bswa, Cswa]
#
# matrizes = []
# pes = []
#
# Perm = eye(3)
#
# def swap_recursi(iter, matriz, perm):
#
#     iter -= 1
#     row = copy.deepcopy(matriz)
#     col = copy.deepcopy(matriz)
#     perm = copy.deepcopy(perm)
#     counter = 0
#
#     for swap in swaps:
#         counter += 1
#         i = swap[0]
#         j = swap[1]
#         row_edit = copy.deepcopy(row)
#         row_edit.row_swap(i,j)
#         perm_edit = copy.deepcopy(perm)
#         perm_edit.row_swap(i,j)
#         matrizes.append(row_edit)
#         pes.append(perm_edit)
#         if iter > 0:
#             swap_recursi(iter,row_edit, perm_edit)
#
#     for swap in swaps:
#         counter += 1
#         i = swap[0]
#         j = swap[1]
#         col_edit = copy.deepcopy(col)
#         col_edit.col_swap(i,j)
#
#         perm_edit = copy.deepcopy(perm)
#         perm_edit.col_swap(i,j)
#
#         matrizes.append(col_edit)
#         pes.append(perm_edit)
#         if iter > 0:
#             swap_recursi(iter,col_edit, perm_edit)
#
#
#
# swap_recursi(4,A,Perm)
#
#
# for matrix in matrizes:
#     a = gauss_seidel(matrix, None, 0, None)
#     Ms.append(a.M)
#
# Ms.sort()
#
# M = min(Ms)
#
#
# print(M)
# last = 100000
# for i in range(0, 500):
#     print(i, M**i/(1-M) * var - erro)

```

```python

```

```python
# x = symbols('x')
# g0 = x*0 + 1
# g1 = x
# g2 = x**2
#
# GS = [g0, g1, g2]
# intervalo = Interval(2,9)
#
# Func = 10/9 + cbrt(2/7*x - 11/7)
#
# a = MMQ_conti(GS,x,Func,intervalo)
#
#

```

max_stars_repo_name: victorathanasio/Personal-projects

max_issues_repo_name: victorathanasio/Personal-projects

max_forks_repo_name: victorathanasio/Personal-projects


text:
Tài liệu này mang giấy phép Creative Commons Attribution (CC BY).
(c) Nguyễn Ngọc Sáng, Zhukovsky 7/2020.

[@SangVn](https://github.com/SangVn) [@VnCFD](https://vncfdgroup.wordpress.com/)

*Thực hành CFD với Python!*

**Danh sách tài liệu tham khảo:**

*    [1. Eleuterio F. Toro] *`Riemann Solvers and Numerical Methods for Fluid Dynamics`*
*    [2. Katake Masatsuka] I *`do like CFD. Governing Equations and Exact Solutions`*
*    [3. F. Moukalled, L. Mangani, M. Darwish] *`The Finite Volume Method in Computational Fluid Dynamics. An Advanced Introduction with OpenFOAM and Matlab`*
*    [4. Randall J. Leveque] *`Finite-Volume Methods for Hyperbolic Problems`*
*    [5. H. K. Versteeg, W. Malalasekera] *`An introduction to Computational Fluid Dynamics. The Finite Volume Method`*

# Bài 25. Hệ phương trình Navier-Stokes, phương pháp Roe tính dòng đối lưu

Hệ phương trình Navier-Stokes được viết rất cụ thể trong quyển sách Abramovich chương 2, đã được dịch và đăng trên trang [VnCFD](https://vncfdgroup.wordpress.com). Các bài mở đầu bao gồm:

* [Bài 1. Chuyển động của phần tử chất lỏng](https://vncfdgroup.wordpress.com/2012/04/09/bai-1-chuy%e1%bb%83n-d%e1%bb%99ng-c%e1%bb%a7a-ph%e1%ba%a7n-t%e1%bb%ad-ch%e1%ba%a5t-l%e1%bb%8fng/)

* [Bài 2. Phương trình liên tục](https://vncfdgroup.wordpress.com/2012/04/14/bai-2-ph%c6%b0%c6%a1ng-trinh-lien-t%e1%bb%a5c/)

* [Bài 3. Các lực tác dụng trong chất lỏng](https://vncfdgroup.wordpress.com/2012/04/14/bai-3-cac-l%e1%bb%b1c-tac-d%e1%bb%a5ng-trong-ch%e1%ba%a5t-l%e1%bb%8fng/)

* [Bài 4. Mối liên hệ giữa ứng suất và tốc độ biến dạng](https://vncfdgroup.wordpress.com/2012/04/14/bai-4-m%e1%bb%91i-lien-h%e1%bb%87-gi%e1%bb%afa-%e1%bb%a9ng-su%e1%ba%a5t-va-d%e1%bb%99-bi%e1%ba%bfn-d%e1%ba%a1ng/)

* [Bài 5. Phương trình Navier – Stokes](https://vncfdgroup.wordpress.com/2012/04/14/bai-5-ph%c6%b0%c6%a1ng-trinh-navier-stokes/)

* [Bài 6. Phương trình năng lượng](https://vncfdgroup.wordpress.com/2012/04/14/bai-6-ph%c6%b0%c6%a1ng-trinh-nang-l%c6%b0%e1%bb%a3ng/)

## 1. Hệ phương trình Navier-Stokes

Để tóm tắt hệ phương trình Navier-Stokes, trước hết ta nhớ lại **hệ phương trình Euler cho khí lý tưởng**:

$$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} + \frac{\partial H}{\partial z} = 0 \qquad (1)$$
với
\begin{align}
U & = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho e\end{pmatrix}, &
F & = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho uv \\ \rho uw \\ \rho uh\end{pmatrix}, &
G & = \begin{pmatrix} \rho v \\ \rho vu \\ \rho v^2 + p \\ \rho vw \\ \rho vh\end{pmatrix}, &
H & = \begin{pmatrix} \rho w \\ \rho wu \\ \rho wv \\ \rho w^2 + p \\ \rho wh\end{pmatrix} &
\end{align}
trong đó: $h = e + \frac{p}{\rho}$, $e = \varepsilon + \frac{u^2+v^2+w^2}{2}$, $\varepsilon = \frac{p}{(\gamma-1)\rho}$; $\rho$ - khối lượng riêng; $u, v, w$ - vận tốc theo phương x, y, z; $p$ - áp suất;  $h$ - enthalpy; $e$ - năng lượng một đơn vị khối lượng chất khí; $\varepsilon$ - nội năng.

Trong phương trình (1), số hạng đầu gọi là **thành phần không ổn định (unsteady term)**, ba số hạng còn lại là **dòng đại lượng (khối lượng, động lượng, năng lượng) đi qua bề mặt thể tích kiểm soát**, được tạo ra bởi chuyển dịch vĩ mô của dòng, gọi là **thành phần đối lưu (convection term)**.

Điểm khác biệt giữa hai hệ phương trình Euler và Navier-Stokes nằm ở chính đối tượng mà chúng mô tả. **Hệ Navier-Stokes mô tả khí thực** do đó phải tính tới lực ma sát trong (lực nhớt), công của lực ma sát, quá trình trao đổi nhiệt (dẫn nhiệt). Hệ này có dạng:

$$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} + \frac{\partial H}{\partial z} = 0 \qquad (2)$$
với
\begin{align}
U & = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho e\end{pmatrix}, &
F & = \begin{pmatrix} \rho u \\ \rho u^2+p\\ \rho uv \\ \rho uw \\ \rho uh\end{pmatrix} +
\begin{pmatrix} 0 \\ \tau_{xx}\\ \tau_{xy} \\ \tau_{xz} \\ u\tau_{xx}+v\tau_{xy}+w\tau_{xz}+q_x\end{pmatrix}, &
\end{align}

\begin{align}
G & = \begin{pmatrix} \rho v \\ \rho vu \\ \rho v^2+p \\ \rho vw \\ \rho vh \end{pmatrix} +
\begin{pmatrix} 0 \\ \tau_{yx} \\ \tau_{yy} \\ \tau_{yz} \\ u\tau_{yx}+v\tau_{yy}+w\tau_{yz}+q_y\end{pmatrix}, &
H & = \begin{pmatrix} \rho w \\ \rho wu \\ \rho wv \\ \rho w^2+p \\ \rho wh \end{pmatrix} +
\begin{pmatrix} 0 \\ \tau_{zx} \\ \tau_{zy} \\ \tau_{zz} \\ u\tau_{zx}+v\tau_{zy}+w\tau_{zz}+q_z\end{pmatrix} &
\end{align}

Trong đó, tensor lực nhớt:
\begin{align}
\tau & = \begin{pmatrix}
\tau_{xx} & \tau_{xy} & \tau_{xz} \\
\tau_{yx} & \tau_{yy} & \tau_{yz} \\
\tau_{zx} & \tau_{zy} & \tau_{zz} \\
\end{pmatrix}
\end{align}
với thành phần theo phương tiếp tuyến:

$$\tau_{xy}=\tau_{yx}=-\mu\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)$$
$$\tau_{yz}=\tau_{zy}=-\mu\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)$$
$$\tau_{zx}=\tau_{xz}=-\mu\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)$$

và thành phần theo phương pháp tuyến:
$$\tau_{xx}=-2\mu\left[\frac{\partial u}{\partial x} - \frac{1}{3}\left(\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)\right]$$

$$\tau_{yy}=-2\mu\left[\frac{\partial v}{\partial y} - \frac{1}{3}\left(\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)\right]$$

$$\tau_{zz}=-2\mu\left[\frac{\partial w}{\partial z} - \frac{1}{3}\left(\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)\right]$$

hệ số nhớt động lực học:
$$\mu = \mu_0\left(\frac{T}{T_0}\right)^{1.5} \frac{T_0+C}{T+C}; T_0=273 [K], C=122 [K], \mu_0=1.72.10^{-5} [Pa.s]. $$

Mật độ dòng nhiệt của quá trình dẫn nhiệt:
$$q_x = -\lambda \frac{\partial T}{\partial x}, q_y = -\lambda \frac{\partial T}{\partial y}, q_z = -\lambda \frac{\partial T}{\partial z},$$
hệ số dẫn nhiệt $\lambda$ tính thông qua số Prandtl:

$$Pr = \frac{c_p \mu}{\lambda}=\frac{\gamma R \mu}{\lambda (\gamma-1)}$$

với không khí trong khoảng nhiệt độ từ 200K tới 1000K, $Pr=0.72$.

Mật độ dòng nhiệt có thể tính theo công thức:

$$q = -\lambda grad T = -\frac{\mu}{Pr(\gamma-1)}grad (c^2) = -\frac{\gamma \mu}{Pr(\gamma-1)}grad (p/\rho), \quad c=\sqrt{\frac{\gamma p}{\rho}}=\sqrt{\gamma R T}$$

So với hệ phương trình Euler, trong Navier-Stokes (2) **dòng qua mặt** gồm hai thành phần: **thành phần đối lưu (convection term)** và **thành phần khuếch tán (diffusion term)**. Thành phần khuếch tán được tạo ra bởi chuyển dịch, va chạm của các phân tử chất khí.

Như vậy, để giải hệ Navier-Stokes ta cần phải thêm vào euler solver **dòng khuếch tán**.

# 2. Phương pháp Roe tính dòng đối lưu

Có nhiều phương pháp tính dòng đối lưu, ở phần 2 *Thực hành CFD với Python* chúng ta đã tìm hiểu phương pháp Godunov: tìm nghiệm biến nguyên thủy của bài toán phân ra gián đoạn Riemann, thông qua đó để tính dòng qua mặt. Ở phần 3, ta đã sử dụng phương pháp Roe để tính dòng qua mặt nhưng chưa tìm hiểu về nó.

Các phương pháp giải bài toán Riemann được viết rất cụ thể trong các quyển sách tham khảo được liệt kê ở trên. Sau đây ta cùng tìm hiểu về phương pháp Roe.

Xét hệ phương trình vi phân:

$$\frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} =0$$

Phương trình này có thể viết lại ở dạng:

$$\frac{\partial U}{\partial t} + A(U)\frac{\partial U}{\partial x} =0$$
với ma trận Jacobi $A(U)=\frac{\partial F}{\partial U}$.

Để giải bài toán Riemann, thay vì sử dụng ma trận chuẩn xác Jacobi A(U) phương pháp xấp xỉ Roe đề xuất sử dụng ma trận hằng số Jacobi $\hat A(U_L, U_R)$. Đối với hệ hyperbolic của $m$ định luật bảo toàn ($m$ phương trình), ma trận $\hat A$ phải thỏa mãn các tính chất sau:

* 1. có các trị riêng là số thực  $\hat\lambda_1 \le \hat\lambda_2 \le... \le \hat\lambda_m$ và
     có hệ đầy đủ vector riêng bên phải độc lập tuyến tính $\hat R_1, \hat R_2, ..., \hat R_m$
* 2. $\hat A(U,U) = A(U)$
* 3. $F(U_R)-F(U_L) = \hat A(U_R-U_L)$

Khi đó ta có:
$$\hat A=\hat R \hat\Lambda \hat R^{-1}$$

với $\hat R, \hat R^{-1}$ là ma trận cột các vector riêng bên phải và nghịch đảo của nó, $\hat \Lambda$ là ma trận đường chéo các trị riêng.

Giá trị của Jacobian $\hat A$ được xác định, sử dụng các giá trị trung bình Roe:

\begin{align}
\hat \rho &= \sqrt{\rho_L \rho_R} \\
\hat u &= \frac{\sqrt{\rho_L}u_L + \sqrt{\rho_R}u_R}{\sqrt{\rho_L}+\sqrt{\rho_R}} \\
\hat v &= \frac{\sqrt{\rho_L}v_L + \sqrt{\rho_R}v_R}{\sqrt{\rho_L}+\sqrt{\rho_R}} \\
\hat h &= \frac{\sqrt{\rho_L}h_L + \sqrt{\rho_R}h_R}{\sqrt{\rho_L}+\sqrt{\rho_R}} \\
\hat c &= \sqrt{(\gamma-1)(\hat h - {\hat q})} \\
\hat q &= (\hat u^2 + \hat v^2)/2
\end{align}

### Công thức tính dòng qua mặt

Trước hết ta biểu diễn hiệu $U_R-U_L$ ở dạng:

$$\Delta U = U_R-U_L = \sum_{i=1}^m \hat \alpha_i \hat R^i$$

Khi đó, giá trị $U_{i+1/2}, F_{i+1/2}$ có thể tính theo 1 trong hai công thức sau:
* Công thức 1:
\begin{align}
U_{i+1/2} &= \frac{1}{2}(U_L + U_R) - \frac{1}{2}\sum_{i=i}^m sign(\hat \lambda_i) \hat \alpha_i \hat R^i \\
F_{i+1/2} &= \frac{1}{2}(F_L + F_R) - \frac{1}{2}\sum_{i=i}^m |\hat \lambda_i| \hat \alpha_i \hat R^i
\end{align}
* Công thức 2:
\begin{align}
U_{i+1/2} &= \frac{1}{2}(U_L + U_R) - \frac{1}{2}R sign(\hat \Lambda) \hat R^{-1} \\
F_{i+1/2} &= \frac{1}{2}(F_L + F_R) - \frac{1}{2}R |\hat \Lambda| \hat R^{-1}
\end{align}

Trong đó $sign(\hat \lambda), |\hat \lambda|$ lần lượt là dấu và trị tuyệt đối của $\hat \lambda$.

Làm thế nào để có công thức tính Jacobian, trị riêng, vector riêng hay $\alpha$?

Xét trường hợp cụ thể sau.
### Bài toán hai chiều
Cụ thể, xét bài toán hai chiều:
\begin{align}
U & = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ \rho e\end{pmatrix}, &
F & = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho uv  \\ \rho uh\end{pmatrix}, &
G & = \begin{pmatrix} \rho v \\ \rho vu \\ \rho v^2 + p \\ \rho vh\end{pmatrix} &
\end{align}

Vector dòng vuông góc với mặt (normal flux) hay là hình chiếu của dòng theo phương pháp tuyến bề mặt $n=(n_x, n_y)^T$ được xác định như sau:

\begin{equation}
F_n = [F,G].n=Fn_x+Gn_y=\begin{pmatrix} \rho V_n \\ \rho V_nu+pn_x \\ \rho V_nv+pn_y \\ \rho V_nh \end{pmatrix}
\end{equation}

với $V_n=un_x+vn_y$. Ma trận Jacobi $A_n=An_x + Bn_y=\frac{\partial F}{\partial U}n_x + \frac{\partial G}{\partial U}n_y$

## 3. Sympy, Jacobian

Để tìm ma trận Jacobi cũng như trị riêng, vector riêng ta sẽ biến đổi toán học bằng sympy.

```python
# Gọi thư viện sympy
from sympy import*

init_printing()
```

```python
# Làm quen với sympy:

# Khai báo các kí tự sẽ sử dụng để biến đổi toán học
r, u, v, p, h, g, q, c = symbols('rho u v p h gamma q c', real=True)
u1, u2, u3, u4 = symbols('u_1 u_2 u_3 u_4', real=True)
nx, ny, Vn = symbols('n_x n_y V_n', real=True)

# Khai báo các vector ở dạng matrận
# q = 0.5(u*u+v*v)

P = Matrix([r, u, v, p])
U = Matrix([r, r*u, r*v, p/(g-1) + r*q])
F = Matrix([r*u, r*u*u + p, r*u*v, (g*p/(g-1) + r*q)*u])
G = Matrix([r*v, r*v*u, r*v*v + p, (g*p/(g-1) + r*q)*v])

# in ra màn hình
P, U, F, G
```

$$\left ( \left[\begin{matrix}\rho\\u\\v\\p\end{matrix}\right], \quad \left[\begin{matrix}\rho\\\rho u\\\rho v\\\frac{p}{\gamma - 1} + q \rho\end{matrix}\right], \quad \left[\begin{matrix}\rho u\\p + \rho u^{2}\\\rho u v\\u \left(\frac{\gamma p}{\gamma - 1} + q \rho\right)\end{matrix}\right], \quad \left[\begin{matrix}\rho v\\\rho u v\\p + \rho v^{2}\\v \left(\frac{\gamma p}{\gamma - 1} + q \rho\right)\end{matrix}\right]\right )$$

Cách xác định ma trận Jacobi (Jacobian):

\begin{equation}
A=\frac{\partial F}{\partial U}=
\frac{\partial(f_1, f_2, f_3, f_4)}{\partial(u_1, u_2, u_3, u_4)}=
\begin{pmatrix}
\frac{\partial f_1}{\partial u_1} & \frac{\partial f_1}{\partial u_2} & \frac{\partial f_1}{\partial u_3} & \frac{\partial f_1}{\partial u_4} \\
\frac{\partial f_2}{\partial u_1} & \frac{\partial f_2}{\partial u_2} & \frac{\partial f_2}{\partial u_3} & \frac{\partial f_2}{\partial u_4} \\
\frac{\partial f_3}{\partial u_1} & \frac{\partial f_3}{\partial u_2} & \frac{\partial f_3}{\partial u_3} & \frac{\partial f_3}{\partial u_4} \\
\frac{\partial f_4}{\partial u_1} & \frac{\partial f_4}{\partial u_2} & \frac{\partial f_4}{\partial u_3} & \frac{\partial f_4}{\partial u_4} \\
\end{pmatrix}
\end{equation}

```python
# Để tìm ma trận Jacobi ta chuyển sang sử dụng các biến (u1,u2,u3,u4)
U = Matrix([u1, u2, u3, u4])

# Khi đó áp suất p được tính theo công thức
pf = (u4-(u2*u2+u3*u3)/(2*u1))*(g-1)

# Biến đổi ta có:
F = Matrix([u2, u2*u2/u1 + pf, u2*u3/u1, (u4+pf)*u2/u1])
G = Matrix([u3, u3*u2/u1, u3*u3/u1 + pf, (u4+pf)*u3/u1])

# Tìm matrận Jacobi
A = F.jacobian(U)
B = G.jacobian(U)

# Để đưa về dạng nguyên bản ta thực hiện bước thay biến
args = ([u1, r], [u2, r*u], [u3, r*v], [u4, p/(g-1)+0.5*r*(u*u+v*v)])

# Dùng hàm simplify để rút gọn công thức, in ra màn hình
A, B = simplify(A.subs(args)), simplify(B.subs(args))
A, B
```

$$\left ( \left[\begin{matrix}0 & 1 & 0 & 0\\- u^{2} + \frac{\left(\gamma - 1\right) \left(u^{2} + v^{2}\right)}{2} & u \left(- \gamma + 3\right) & v \left(- \gamma + 1\right) & \gamma - 1\\- u v & v & u & 0\\\frac{u \left(\gamma^{2} \rho u^{2} + \gamma^{2} \rho v^{2} - 2 \gamma p - 3 \gamma \rho u^{2} - 3 \gamma \rho v^{2} + 2 \rho u^{2} + 2 \rho v^{2}\right)}{2 \rho \left(\gamma - 1\right)} & \frac{- 2 \gamma^{2} \rho u^{2} + 2 \gamma p + 5 \gamma \rho u^{2} + \gamma \rho v^{2} - 3 \rho u^{2} - \rho v^{2}}{2 \rho \left(\gamma - 1\right)} & u v \left(- \gamma + 1\right) & \gamma u\end{matrix}\right], \quad \left[\begin{matrix}0 & 0 & 1 & 0\\- u v & v & u & 0\\- v^{2} + \frac{\left(\gamma - 1\right) \left(u^{2} + v^{2}\right)}{2} & u \left(- \gamma + 1\right) & v \left(- \gamma + 3\right) & \gamma - 1\\\frac{v \left(\gamma^{2} \rho u^{2} + \gamma^{2} \rho v^{2} - 2 \gamma p - 3 \gamma \rho u^{2} - 3 \gamma \rho v^{2} + 2 \rho u^{2} + 2 \rho v^{2}\right)}{2 \rho \left(\gamma - 1\right)} & u v \left(- \gamma + 1\right) & \frac{- 2 \gamma^{2} \rho v^{2} + 2 \gamma p + \gamma \rho u^{2} + 5 \gamma \rho v^{2} - \rho u^{2} - 3 \rho v^{2}}{2 \rho \left(\gamma - 1\right)} & \gamma v\end{matrix}\right]\right )$$

```python
# Jacobian An
An = A*nx+B*ny

# Tìm trị riêng
simplify(An.eigenvals())
```

Như vậy Jacobian An có 4 trị riêng:
$$\lambda_1=\lambda_2=V_n, \lambda_3=V_n-c, \lambda_4=V_n+c$$
với vận tốc âm thanh $c=\sqrt{\gamma p/ \rho}$

Để tìm vector riêng bên phải, bên trái có thể dùng hàm **An.eigenvects(), An.left_eigenvects()**, tuy nhiên công thức thu được sẽ rất rất dài vì không được rút gọn. Vì vậy ta sẽ thu gọn Jacobian An và tính toán thủ công.

```python
# Rút gọn công thức ta có hai ma trận Jacobi sau:
A = Matrix([[0, 1, 0, 0], [(g-1)*q-u*u, (3-g)*u, (1-g)*v, g-1],\
            [-u*v, v, u, 0], [((g-1)*q-h)*u, h+(1-g)*u*u, (1-g)*u*v, g*u]])

B = Matrix([[0, 0, 1, 0], [-u*v, v, u, 0], \
                [(g-1)*q-v*v, (1-g)*u, (3-g)*v, g-1],\
                [((g-1)*q-h)*v, (1-g)*u*v, h+(1-g)*v*v, g*v]])
# In ra màn hình
# A, B
```

```python
# Jacobian An
An = A*nx+B*ny
An
```

$$\left[\begin{matrix}0 & n_{x} & n_{y} & 0\\n_{x} \left(q \left(\gamma - 1\right) - u^{2}\right) - n_{y} u v & n_{x} u \left(- \gamma + 3\right) + n_{y} v & n_{x} v \left(- \gamma + 1\right) + n_{y} u & n_{x} \left(\gamma - 1\right)\\- n_{x} u v + n_{y} \left(q \left(\gamma - 1\right) - v^{2}\right) & n_{x} v + n_{y} u \left(- \gamma + 1\right) & n_{x} u + n_{y} v \left(- \gamma + 3\right) & n_{y} \left(\gamma - 1\right)\\n_{x} u \left(- h + q \left(\gamma - 1\right)\right) + n_{y} v \left(- h + q \left(\gamma - 1\right)\right) & n_{x} \left(h + u^{2} \left(- \gamma + 1\right)\right) + n_{y} u v \left(- \gamma + 1\right) & n_{x} u v \left(- \gamma + 1\right) + n_{y} \left(h + v^{2} \left(- \gamma + 1\right)\right) & \gamma n_{x} u + \gamma n_{y} v\end{matrix}\right]$$

```python
# Kiểm tra lại trị riêng
An.eigenvals()
```

Rút gọn ta được 4 trị riêng như ở trên:
$$\lambda_1=\lambda_2=V_n, \lambda_3=V_n-c, \lambda_4=V_n+c$$

```python
# Rút gọn công thức ta có Jacobian An:
An = Matrix([[0 , nx , ny , 0], \
            [(g-1)*q*nx - Vn*u, (2-g)*u*nx + Vn, (1-g)*v*nx + u*ny , (g-1)*nx], \
            [(g-1)*q*ny - Vn*v, (1-g)*u*ny + v*nx, (2-g)*v*ny + Vn , (g-1)*ny], \
            [((g-1)*q - h)*Vn, (1-g)*Vn*u + h*nx , (1-g)*Vn*v + h*ny, g*Vn]])

# Tìm vector riêng bên phải:
# 1. trường hợp trị riêng bằng Vn
A1 = An - Vn*eye(4) # = An - lambda*I
A1
```

$$\left[\begin{matrix}- V_{n} & n_{x} & n_{y} & 0\\- V_{n} u + n_{x} q \left(\gamma - 1\right) & n_{x} u \left(- \gamma + 2\right) & n_{x} v \left(- \gamma + 1\right) + n_{y} u & n_{x} \left(\gamma - 1\right)\\- V_{n} v + n_{y} q \left(\gamma - 1\right) & n_{x} v + n_{y} u \left(- \gamma + 1\right) & n_{y} v \left(- \gamma + 2\right) & n_{y} \left(\gamma - 1\right)\\V_{n} \left(- h + q \left(\gamma - 1\right)\right) & V_{n} u \left(- \gamma + 1\right) + h n_{x} & V_{n} v \left(- \gamma + 1\right) + h n_{y} & V_{n} \gamma - V_{n}\end{matrix}\right]$$

Vì trị riêng Vn lặp lại hai lần nên tương ứng ta lấy hai vector riêng bên phải. Dễ dàng dàng nhận thấy hai vector này có dạng:

$R_1 = (1, u, v, \tau_1)^T, R_2 = (0, -ny, nx, \tau_1)^T$

Chúng ta cần tìm $\tau_1$ cho từng trường hợp. Để làm điều này ta sử dụng hàm gauss_jordan_solve để hỗ trợ. Cách làm như sau:

```python
B = Matrix([0, 0, 0, 0])

Ri, params = A1.gauss_jordan_solve(B)
simplify(Ri)
```

$$\left[\begin{matrix}\frac{n_{x} \tau_{0} v - n_{x} \tau_{1} - n_{y} \tau_{0} u}{- V_{n} u + n_{x} q}\\\frac{V_{n} \tau_{0} v - V_{n} \tau_{1} - n_{y} q \tau_{0}}{- V_{n} u + n_{x} q}\\\tau_{0}\\\tau_{1}\end{matrix}\right]$$

Đồng nhất hạng tử Ri=R1, Ri=R2 ta tìm được $\tau_1$ lần lượt bằng $q; n_xv-n_yu$. Như vậy hai vector riêng bên phải:

$$R_1 = (1, u, v, q)^T$$ $$R_2 = (0, -ny, nx, n_xv-n_yu)^T$$

```python
# 2. trường hợp trị riêng  bằng Vn-c
A2 = An - (Vn-c)*eye(4)
A2
```

$$\left[\begin{matrix}- V_{n} + c & n_{x} & n_{y} & 0\\- V_{n} u + n_{x} q \left(\gamma - 1\right) & c + n_{x} u \left(- \gamma + 2\right) & n_{x} v \left(- \gamma + 1\right) + n_{y} u & n_{x} \left(\gamma - 1\right)\\- V_{n} v + n_{y} q \left(\gamma - 1\right) & n_{x} v + n_{y} u \left(- \gamma + 1\right) & c + n_{y} v \left(- \gamma + 2\right) & n_{y} \left(\gamma - 1\right)\\V_{n} \left(- h + q \left(\gamma - 1\right)\right) & V_{n} u \left(- \gamma + 1\right) + h n_{x} & V_{n} v \left(- \gamma + 1\right) + h n_{y} & V_{n} \gamma - V_{n} + c\end{matrix}\right]$$

```python
# Không còn dùng được hàm gauss_jordan_solve để giải
A2.gauss_jordan_solve(B)
```

$$\left ( \left[\begin{matrix}0\\0\\0\\0\end{matrix}\right], \quad \left[\begin{matrix}\end{matrix}\right]\right )$$

```python
# Vector riêng có dạng Ri = (1, u-c*nx, v-c*ny, x).
x = symbols('x', real=True)
Ri = Matrix([1, u-c*nx, v-c*ny, x])

# Tìm x như sau:
fx = A2[3, :]*Ri # = 0

solve(fx, x)
```

Rút gọn tử số ta được: $-V_n^2c(\gamma-1)+V_nq(\gamma-1)+ch$. Thay $q=h-c^2/(\gamma-1)$, tử số bằng:

$-V_n^2c(\gamma-1) + V_n(\gamma-1)h - V_nc^2+ch = (V_n(\gamma-1)+c)(h-V_nc)$

Như vậy $x=h-V_nc$.

Ta có vector riêng tương ứng trị riêng $\lambda = V_n - c$:
$$R_3=(1, u-cn_x, v-cn_y, h-V_nc)^T$$
Tương tự, vector riêng tương ứng trị riêng $\lambda = V_n + c$:
$$R_4 = (1, u+cn_x, v+cn_y, n+V_nc)^T$$

```python
# Kiểm tra lại vector riêng R3
Ri = Matrix([1, u-c*nx, v-c*ny, h-Vn*c])
fx = A2[2, :]*Ri

args = ([Vn, u*nx+v*ny], [q, 0.5*(u*u+v*v)], [h, 0.5*(u*u+v*v) + c*c/(g-1)])
simplify(fx.subs(args)) # = 0 -> đúng
```

$$\left[\begin{matrix}c v \left(- n_{x}^{2} - n_{y}^{2} + 1\right)\end{matrix}\right]$$

```python
fx = A2[1, :]*Ri
simplify(fx.subs(args)) # = 0 -> đúng
```

$$\left[\begin{matrix}c u \left(- n_{x}^{2} - n_{y}^{2} + 1\right)\end{matrix}\right]$$

```python
# Ta có ma trận đường chéo trị riêng
E = diag(abs(Vn-c), abs(Vn), abs(Vn), abs(Vn+c))
E
```

$$\left[\begin{matrix}\left|{V_{n} - c}\right| & 0 & 0 & 0\\0 & \left|{V_{n}}\right| & 0 & 0\\0 & 0 & \left|{V_{n}}\right| & 0\\0 & 0 & 0 & \left|{V_{n} + c}\right|\end{matrix}\right]$$

```python
# Tương ứng ma trận vector riêng bên phải
R = Matrix([[1, 1, 0, 1], [u-c*nx, u, -ny, u+c*nx],\
            [v-c*ny, v, nx, v+c*ny], [h-Vn*c, q, -u*ny+v*nx, h+Vn*c]])

R
```

$$\left[\begin{matrix}1 & 1 & 0 & 1\\- c n_{x} + u & u & - n_{y} & c n_{x} + u\\- c n_{y} + v & v & n_{x} & c n_{y} + v\\- V_{n} c + h & q & n_{x} v - n_{y} u & V_{n} c + h\end{matrix}\right]$$

Kết quả này hoàn toàn trùng khớp với kết quả trong sách của **Katate Masatsuka** trang 69. Tương tự ta có ma trận vector riêng bên trái:

```python
L = Matrix([[((g-1)*q/c**2 + Vn/c)/2, -((g-1)*u/c**2 + nx/c)/2, -((g-1)*v/c**2+ny/c)/2, (g-1)/(2*c**2)], \
            [1-(g-1)*q/c**2, (g-1)*u/c**2, (g-1)*v/c**2, -(g-1)/c**2], \
            [u*ny-v*nx, -ny, nx, 0], \
            [((g-1)*q/c**2-Vn/c)/2, -((g-1)*u/c**2-nx/c)/2, -((g-1)*v/c**2-ny/c)/2, (g-1)/(2*c**2)]])

L
```

$$\left[\begin{matrix}\frac{V_{n}}{2 c} + \frac{q \left(\gamma - 1\right)}{2 c^{2}} & - \frac{n_{x}}{2 c} - \frac{u \left(\gamma - 1\right)}{2 c^{2}} & - \frac{n_{y}}{2 c} - \frac{v \left(\gamma - 1\right)}{2 c^{2}} & \frac{\gamma - 1}{2 c^{2}}\\1 - \frac{q \left(\gamma - 1\right)}{c^{2}} & \frac{u \left(\gamma - 1\right)}{c^{2}} & \frac{v \left(\gamma - 1\right)}{c^{2}} & \frac{- \gamma + 1}{c^{2}}\\- n_{x} v + n_{y} u & - n_{y} & n_{x} & 0\\- \frac{V_{n}}{2 c} + \frac{q \left(\gamma - 1\right)}{2 c^{2}} & \frac{n_{x}}{2 c} - \frac{u \left(\gamma - 1\right)}{2 c^{2}} & \frac{n_{y}}{2 c} - \frac{v \left(\gamma - 1\right)}{2 c^{2}} & \frac{\gamma - 1}{2 c^{2}}\end{matrix}\right]$$

```python
# Kiểm tra lại kết quả, tính tích R*L
RL = simplify(R*L)
simplify(RL.subs(args))
```

$$\left[\begin{matrix}1 & 0 & 0 & 0\\u \left(- n_{x}^{2} - n_{y}^{2} + 1\right) & n_{x}^{2} + n_{y}^{2} & 0 & 0\\v \left(- n_{x}^{2} - n_{y}^{2} + 1\right) & 0 & n_{x}^{2} + n_{y}^{2} & 0\\- 1.0 n_{x}^{2} u^{2} - 1.0 n_{x}^{2} v^{2} - 1.0 n_{y}^{2} u^{2} - 1.0 n_{y}^{2} v^{2} + 1.0 u^{2} + 1.0 v^{2} & u \left(n_{x}^{2} + n_{y}^{2} - 1\right) & v \left(n_{x}^{2} + n_{y}^{2} - 1\right) & 1\end{matrix}\right]$$

Rút gọn ta được $RL=I$, hay $L=R^{-1}$.

Như vậy ta đã có công thức của $R, \Lambda, L$, để tính dòng theo công thức 1, cần xác định công thức tính $\alpha_i$. Cách làm như trong sách của **Toro, trang 364**:

\begin{align}
\hat\alpha_1 &= \frac{1}{2\hat c^2}[\Delta p-\hat \rho \hat c \Delta V_n]\\
\hat\alpha_2 &= \Delta \rho -\Delta p/\hat c^2 \\
\hat\alpha_3 &= \hat \rho \Delta V_{\tau}\\
\hat\alpha_4 &= \frac{1}{2\hat c^2}[\Delta p+\hat \rho \hat c \Delta V_n]\\
\end{align}

với $V_n = Vn = (u,v)(n_x,n_y) = un_x+vn_y; \quad V_{\tau} = V{\tau} = (u,v)(-n_y,n_x) = -un_y+vn_x$.

## 5. So sánh hai công thức tính dòng qua mặt

Viết hai chương trình nhỏ để kiểm tra phương pháp tính dòng.

```python
from numpy import diag, array, dot
```

```python
# Công thức 1
def flux_roe_1(Pl, Pr, N):
    # Xác định các thông số trái, phải
    rl, ul, vl, pl = Pl[0], Pl[1], Pl[2], Pl[3]
    rr, ur, vr, pr = Pr[0], Pr[1], Pr[2], Pr[3]
    nx, ny = N[0], N[1]

    g, gm1 = 1.4, 0.4
    cl = (g*pl/rl)**0.5
    cr = (g*pr/rr)**0.5
    hl = 0.5*(ul*ul+vl*vl)+cl*cl/gm1
    hr = 0.5*(ur*ur+vr*vr)+cr*cr/gm1
    Vln = ul*nx+vl*ny
    Vlt = -ul*ny+vl*nx
    Vrn = ur*nx+vr*ny
    Vrt = -ur*ny+vr*nx

    # Xác định các giá trị trung bình Roe
    rla = rl**0.5
    rra = rr**0.5
    rsa = rla+rra

    r = rla*rra
    u = (rla*ul+rra*ur)/rsa
    v = (rla*vl+rra*vr)/rsa
    Vn = u*nx+v*ny
    h = (rla*hl+rra*hr)/rsa
    q = 0.5*(u*u+v*v)
    c = (gm1*(h-q))**0.5
    c2 = c*c

    # Độ lớn Delta U (để tính alpha)
    Dr  = rr-rl
    DVn = Vrn-Vln
    DVt = Vrt-Vlt
    Dp  = pr-pl

    # Ma trận vector riêng bên phải
    R = array([[1, 1, 0, 1], [u-c*nx, u, -ny, u+c*nx],\
               [v-c*ny, v, nx, v+c*ny], [h-Vn*c, q, -u*ny+v*nx, h+Vn*c]])

    # Vận tốc sóng
    E = [abs(Vn-c), abs(Vn), abs(Vn), abs(Vn+c)]

    # Alpha
    A = [(Dp-r*c*DVn)/(2*c2), Dr-Dp/c2, r*DVt, (Dp+r*c*DVn)/(2*c2)]

    dF = array([0., 0., 0., 0.])
    for i in range(4):
        for j in range(4):
            dF[i] = dF[i] + E[j]*A[j]*R[i,j]

    # Dòng trái, phải
    rln = rl*Vln
    Fl  = array([rln, rln*ul+pl*nx, rln*vl+pl*ny, rln*hl])
    rrn = rr*Vrn
    Fr  = array([rrn, rrn*ur+pr*nx, rrn*vr+pr*ny, rrn*hr])

    # Dòng qua mặt
    F = 0.5*(Fl+Fr - dF)

    return F
```

```python
# Công thức 2
def flux_roe_2(Pl, Pr, N):
    # Xác định các thông số trái, phải
    rl, ul, vl, pl = Pl[0], Pl[1], Pl[2], Pl[3]
    rr, ur, vr, pr = Pr[0], Pr[1], Pr[2], Pr[3]
    nx, ny = N[0], N[1]

    g, gm1 = 1.4, 0.4
    cl = (g*pl/rl)**0.5
    cr = (g*pr/rr)**0.5
    hl = 0.5*(ul*ul+vl*vl)+cl*cl/gm1
    hr = 0.5*(ur*ur+vr*vr)+cr*cr/gm1
    Vln = ul*nx+vl*ny
    Vrn = ur*nx+vr*ny

    # Xác định các giá trị trung bình Roe
    rla = rl**0.5
    rra = rr**0.5
    rsa = rla+rra

    r = rla*rra
    u = (rla*ul+rra*ur)/rsa
    v = (rla*vl+rra*vr)/rsa
    Vn = u*nx+v*ny
    h = (rla*hl+rra*hr)/rsa
    q = 0.5*(u*u+v*v)
    c = (gm1*(h-q))**0.5
    gm1dc2 = gm1/(c*c)

    # Ma trận vector riêng bên phải
    R = array([[1, 1, 0, 1], [u-c*nx, u, -ny, u+c*nx],\
               [v-c*ny, v, nx, v+c*ny], [h-Vn*c, q, -u*ny+v*nx, h+Vn*c]])

    # Ma trận đường chéo trị riêng
    E = diag([abs(Vn-c), abs(Vn), abs(Vn), abs(Vn+c)])

    # Ma trận vector riêng bên trái
    L = array([[(gm1dc2*q + Vn/c)/2, -(gm1dc2*u + nx/c)/2, -(gm1dc2*v+ny/c)/2, gm1dc2/2], \
               [1-gm1dc2*q, gm1dc2*u, gm1dc2*v, -gm1dc2], \
               [u*ny-v*nx, -ny, nx, 0],\
               [(gm1dc2*q-Vn/c)/2, -(gm1dc2*u-nx/c)/2, -(gm1dc2*v-ny/c)/2, gm1dc2/2]])

    # Tích R*E*L
    RE = dot(R,E)
    REL = dot(RE,L)

    # Biến bảo toàn trái, phải
    Ul = array([rl, rl*ul, rl*vl, rl*(hl-pl/rl)])
    Ur = array([rr, rr*ur, rr*vr, rr*(hr-pr/rr)])

    # Dòng trái, phải
    rln = rl*Vln
    Fl = array([rln, rln*ul+pl*nx, rln*vl+pl*ny, rln*hl])
    rrn = rr*Vrn
    Fr = array([rrn, rrn*ur+pr*nx, rrn*vr+pr*ny, rrn*hr])

    # Dòng qua mặt
    F = 0.5*(Fl+Fr - dot(REL,Ur-Ul))

    return F
```

```python
# Hàm test thời gian tính toán
def test(flux_func):
    for i in range(10000): F=flux_func(Pl,Pr,N)
    print(F)
    return
```

```python
# Ví dụ
Pl = [1.0, 0.75, 0.0, 1.0]
Pr = [0.125, 0.0, 0.0, 0.1]
N  = [1, 0]
```

```python
%time test(flux_roe_1)
```

    [0.88328704 1.4815703  0.         3.22000163]
    CPU times: user 197 ms, sys: 0 ns, total: 197 ms
    Wall time: 196 ms

```python
%time test(flux_roe_2)
```

    [0.88328704 1.4815703  0.         3.22000163]
    CPU times: user 243 ms, sys: 0 ns, total: 243 ms
    Wall time: 243 ms

Kết quả cho thấy có thể sử dụng cả hai công thức để tính dòng qua mặt, tuy nhiên công thức thứ 2 sẽ chạy lâu hơn.

Còn nhiều phương pháp khác để tính dòng qua mặt được giới thiệu trong sách của Toro. Các bạn có thể tham khảo thêm.

## 6. Entropy fix

Một vấn đề có thể gặp phải khi sử dụng phương pháp Roe đó là kết quả thu được có thể không thỏa mãn điều kiện entropy, lý do là vì nghiệm thu được của bài toán Riemann bằng phương pháp Roe chỉ chứa gián đoạn tiếp xúc hoặc sóng xung kích mà không có sóng giãn **[Toro, trang 366]**. Có nhiều cách để giải quyết vấn đề này, một trong số đó là **Harten's Entropy Fix**: thay vì sử dụng giá trị tuyệt đối $|\lambda|$ ta dùng hàm
$$
\phi_{\delta}(\lambda) =
\begin{cases}
|\lambda| \ \qquad \qquad if \ |\lambda| \ge \delta \\
(\lambda^2+\delta^2)/(2\delta) \ if \ |\lambda| < \delta
\end{cases}
$$
với, ví dụ $\delta=0.2$.

Để ý rằng, sóng giãn chỉ có thể xuất hiện với hai trị riêng $V_n-c, V_n+c$ (hai giá trị còn lại tương ứng với gián đoạn tiếp xúc) do đó ta chỉ cần áp dụng entropy fix với hai trị riêng này.

# 7. Tổng kết

Như vậy, ta đã tìm hiểu được hai phương pháp tính dòng đối lưu qua mặt, đó là phương pháp Godunov và phương pháp Roe. Còn rất nhiều phương pháp khác, các bạn hãy tìm hiểu thêm. Ở bài tiếp theo sẽ giới thiệu về phương pháp tính dòng khuếch tán qua mặt.

# [Bài 26. Phương pháp tính dòng khuếch tán ](Bai_26.ipynb)

```python

```


text:
## 4.6 概率论与数理统计

概率论与数理统计是许多方向开展科学研究的基础，不管是描述客观存在的数据，还是刻画变量之间的关联规则，凭借概率论知识都能得心应手。概率论的数学公式也有其自身特点，本节主要介绍概率论与数理统计范畴内常用的数学公式在LaTeX中的写法。

### 4.6.1 概率论基础

概率论中有一个重要的准则叫做贝叶斯准则，贝叶斯准则的基础为贝叶斯公式，用来描述两个条件概率之间的关系。

【**例1**】书写贝叶斯公式$p\left(\theta\mid y\right)
=\frac{p\left(\theta,y\right)}{p\left(y\right)}
=\frac{p\left(\theta\right)p\left(y\mid\theta\right)}{p\left(y\right)}$和$p\left(\theta\mid y\right)\propto p\left(\theta\right)p\left(y\mid\theta\right)$。

```tex
\documentclass[12pt]{article}
\begin{document}

$$p\left(\theta\mid y\right)
=\frac{p\left(\theta,y\right)}{p\left(y\right)}
=\frac{p\left(\theta\right)p\left(y\mid\theta\right)}{p\left(y\right)}$$
$$p\left(\theta\mid y\right)\propto p\left(\theta\right)p\left(y\mid\theta\right)$$

\end{document}
```
根据贝叶斯公式中的思想，许多机器学习模型的参数可以用贝叶斯算法来求解。

【**例2**】书写$p\left(y\right)=\int p\left(y,\theta\right)d\theta=\int p\left(\theta\right)p\left(y\mid\theta\right)d\theta$。

```tex
\documentclass[12pt]{article}
\begin{document}

$$p\left(y\right)=\int p\left(y,\theta\right)\,d\theta
=\int p\left(\theta\right)p\left(y\mid\theta\right)\,d\theta$$

\end{document}
```

【**例3**】书写期望$\mathbb{E}\left(x\right)=\int xp\left(x\right)dx$与方差$\mathbb{V}\left(x\right)=\int\left(x-\mathbb{E}\left(x\right)\right)^{2}p\left(x\right)dx$。

```tex
\documentclass[12pt]{article}
\begin{document}

$$\mathbb{E}\left(x\right)=\int xp\left(x\right)dx$$
$$\mathbb{V}\left(x\right)=
\int\left(x-\mathbb{E}\left(x\right)\right)^{2}p\left(x\right)dx$$

\end{document}
```

### 4.6.2 概率分布

概率分布，是指用于表述随机变量取值的概率规律，是概率论中最常见的表达式之一。

【**例4**】书写正态分布$x\sim\mathcal{N}\left(\mu,\sigma^{2}\right)$及$p\left(x\right)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2\sigma^{2}}\left(x-\mu\right)^{2}\right)$。

```tex
\documentclass[12pt]{article}
\begin{document}

$$x\sim\mathcal{N}\left(\mu,\sigma^{2}\right)$$
$$p\left(x\right)=\frac{1}{\sqrt{2\pi}\sigma}
\exp\left(-\frac{1}{2\sigma^{2}}\left(x-\mu\right)^{2}\right)$$

\end{document}
```

【**例5**】书写

\begin{equation}
p(y)=\frac{\operatorname{Poisson}(y\mid\theta)
\operatorname{Gamma}(\theta\mid\alpha,\beta)}
{\operatorname{Gamma}(\theta\mid\alpha+y,1+\beta)}
=\frac{\Gamma(\alpha+y)\beta^{\alpha}}
{\Gamma(\alpha)y!(1+\beta)^{\alpha+y}}
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

\begin{equation}
p(y)=\frac{\operatorname{Poisson}(y\mid\theta)
\operatorname{Gamma}(\theta\mid\alpha,\beta)}
{\operatorname{Gamma}(\theta\mid\alpha+y,1+\beta)}
=\frac{\Gamma(\alpha+y)\beta^{\alpha}}
{\Gamma(\alpha)y!(1+\beta)^{\alpha+y}}
\end{equation}

\end{document}
```

【**例6**】书写

$$\theta\mid y\sim\operatorname{Gamma}\left(\alpha+\sum_{i=1}^{n}y_{i},\beta+\sum_{i=1}^{n}x_{i}\right)$$

```tex
\documentclass[12pt]{article}
\begin{document}

$$\theta\mid y\sim\operatorname{Gamma}\left(\alpha+\sum_{i=1}^{n}y_{i},\beta+\sum_{i=1}^{n}x_{i}\right)$$

\end{document}
```

【**例7**】书写

$$\sigma^{2}\mid y\sim\operatorname{Inv}-\chi^{2}\left(n,s^{2}\right)$$

```tex
\documentclass[12pt]{article}
\begin{document}

$$\sigma^{2}\mid y\sim\operatorname{Inv}-\Chi^{2}\left(n,s^{2}\right)$$

\end{document}
```

【**例8**】书写

$$p\left(\beta\mid\mu_{1},\mu_{2},\tau_{1},\tau_{2},\rho\right)
=\prod_{j=1}^{J}\mathcal{N}\left(\begin{pmatrix}\beta_{1j} \\ \beta_{2j} \end{pmatrix}
\bigg|\begin{pmatrix} \mu_{1} \\ \mu_{2} \end{pmatrix},
\begin{pmatrix} \tau_{1}^{2} & \rho\tau_{1}\tau_{2} \\
\rho\tau_{1}\tau_{2} & \tau_{2}^{2} \end{pmatrix}\right)$$

```tex
\documentclass[12pt]{article}
\begin{document}

$$p\left(\beta\mid\mu_{1},\mu_{2},\tau_{1},\tau_{2},\rho\right)
=\prod_{j=1}^{J}\mathcal{N}\left(\begin{pmatrix}\beta_{1j} \\ \beta_{2j} \end{pmatrix}
\bigg|\begin{pmatrix} \mu_{1} \\ \mu_{2} \end{pmatrix},
\begin{pmatrix} \tau_{1}^{2} & \rho\tau_{1}\tau_{2} \\
\rho\tau_{1}\tau_{2} & \tau_{2}^{2} \end{pmatrix}\right)$$

\end{document}
```

【**例9**】书写

\begin{equation}
\begin{aligned}
y_{ij}&\sim\mathcal{N}\left(\alpha_{j}+x_{ij}\beta_{j},\sigma_{y}^{2}\right) \\
\begin{pmatrix}\alpha \\ \beta \end{pmatrix}&
\sim\mathcal{N}\left(\begin{pmatrix} \mu_{\alpha} \\
\mu_{\beta} \end{pmatrix},
\begin{pmatrix} \sigma_{\alpha}^{2} &
\rho\sigma_{\alpha}\sigma_{\beta} \\
\rho\sigma_{\alpha}\sigma_{\beta} &
\sigma_{\beta}^{2} \end{pmatrix}\right)
\end{aligned}
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

\begin{equation}
\begin{aligned}
y_{ij}&\sim\mathcal{N}\left(\alpha_{j}+x_{ij}\beta_{j},\sigma_{y}^{2}\right) \\
\begin{pmatrix}\alpha \\ \beta \end{pmatrix}&
\sim\mathcal{N}\left(\begin{pmatrix} \mu_{\alpha} \\
\mu_{\beta} \end{pmatrix},
\begin{pmatrix} \sigma_{\alpha}^{2} &
\rho\sigma_{\alpha}\sigma_{\beta} \\
\rho\sigma_{\alpha}\sigma_{\beta} &
\sigma_{\beta}^{2} \end{pmatrix}\right)
\end{aligned}
\end{equation}

\end{document}
```

### 练习题

> 打开LaTeX在线系统[https://www.overleaf.com](https://www.overleaf.com/project)或本地安装好的LaTeX编辑器，创建名为LaTeX_practice的项目，并同时新建一个以`.tex`为拓展名的源文件，完成以下几道练习题。

1. 书写

\begin{equation}
p\left(\alpha,\beta\mid y\right)\propto p\left(\alpha,\beta\right)\prod_{j=1}^{J}\frac{\Gamma\left(\alpha+\beta\right)}{\Gamma\left(\alpha\right)\Gamma\left(\beta\right)}\frac{\Gamma\left(\alpha+y_{j}\right)\Gamma\left(\beta+n_{j}-y_{j}\right)}{\Gamma\left(\alpha+\beta+n_{j}\right)}
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

%% 提示：Gamma函数对应的命令为\Gamma
\begin{equation}
% 在此处书写公式
\end{equation}

\end{document}
```

2. 书写

\begin{equation}
\text{Pr}\left(y_{jt}=1\right)=\operatorname{logit}^{-1}\left(\beta_{0}+\beta_{1}x_{j1}+\beta_{2}x_{j2}+\beta_{3}\left(1-x_{j2}\right)t+\beta_{4}x_{j2}t+\alpha_{j}\right)
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

%% 提示：logit函数对应的命令为\operatorname{logit}
\begin{equation}
% 在此处书写公式
\end{equation}

\end{document}
```

3. 书写

\begin{equation}
p_{i}=\text{Pr}\left(y_{i}^{\text{rep}}<y_{i}\mid y\right)+\frac{1}{2}\text{Pr}\left(y_{i}^{\text{rep}}=y_{i}\mid y\right)
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

\begin{equation}
% 在此处书写公式
\end{equation}

\end{document}
```

4. 书写

\begin{equation}
p\left(\theta\mid x,y_{\text{obs}},I\right)=p\left(\theta\mid x\right)\iint p\left(\phi\mid x,\theta\right)p\left(y\mid x,\theta\right)p\left(I\mid x,y,\phi\right)\,\mathrm{d}y_{\text{mis}}\mathrm{d}\phi
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

\begin{equation}
% 在此处书写公式
\end{equation}

\end{document}
```

5. 书写

\begin{equation}
\log p\left(\theta\mid y\right)=-\frac{1}{2}\sum_{j=1}^{J}\frac{\left(y_{j}-\alpha_{j}\right)^{2}}{\sigma_{j}^{2}}-8\log\tau-\frac{1}{2\tau^{2}}\sum_{j=1}^{J}\left(\alpha_{j}-\mu\right)^{2}+\text{const}
\end{equation}

```tex
\documentclass[12pt]{article}
\begin{document}

%% 提示：log函数对应的命令为\log
\begin{equation}
% 在此处书写公式
\end{equation}

\end{document}
```

【回放】[**4.5 线性代数**](https://nbviewer.jupyter.org/github/xinychen/latex-cookbook/blob/main/chapter-4/section5.ipynb)

【继续】[**4.7 优化理论**](https://nbviewer.jupyter.org/github/xinychen/latex-cookbook/blob/main/chapter-4/section7.ipynb)

### License

This work is released under the MIT license.

max_stars_repo_name: Enidsky/latex-cookbook

max_issues_repo_name: liudmmmmm/latex-cookbook

max_forks_repo_name: liudmmmmm/latex-cookbook


text:
## Pattern Recognition
### Assignment 4
#### Group 4:

- COE18B056 - Thigulla Vamsi Krishna
- COE18B065 - Srinivasan R Sharma
- CED18I039 - Paleti Krishnasai

**Q4.)
From the iris dataset, choose the ’petallength’, ’sepalwidth’ for setosa, versicolor and virginica flowers.  Learna decision boundary for the two features using asingle perceptron and SVM. Assume that all the weightsof the perceptron are initialized as 0 with the learning rate of 0.01.  Draw the decision boundary**

```python
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import math
from cvxopt import matrix, solvers
from sympy import *
```

```python
iris_dataset = pd.read_csv("Iris_dataset.csv",index_col=False)
print(iris_dataset.head())
```

       sepal.length  sepal.width  petal.length  petal.width variety
    0           5.1          3.5           1.4          0.2  Setosa
    1           4.9          3.0           1.4          0.2  Setosa
    2           4.7          3.2           1.3          0.2  Setosa
    3           4.6          3.1           1.5          0.2  Setosa
    4           5.0          3.6           1.4          0.2  Setosa

```python
iris_data = iris_dataset.drop(['sepal.length','petal.width'], axis=1)
print(iris_data)
```

         sepal.width  petal.length    variety
    0            3.5           1.4     Setosa
    1            3.0           1.4     Setosa
    2            3.2           1.3     Setosa
    3            3.1           1.5     Setosa
    4            3.6           1.4     Setosa
    ..           ...           ...        ...
    145          3.0           5.2  Virginica
    146          2.5           5.0  Virginica
    147          3.0           5.2  Virginica
    148          3.4           5.4  Virginica
    149          3.0           5.1  Virginica

    [150 rows x 3 columns]

```python
iris_data_12 = iris_data[iris_data['variety']!='Virginica']
print(iris_data_12)
```

        sepal.width  petal.length     variety
    0           3.5           1.4      Setosa
    1           3.0           1.4      Setosa
    2           3.2           1.3      Setosa
    3           3.1           1.5      Setosa
    4           3.6           1.4      Setosa
    ..          ...           ...         ...
    95          3.0           4.2  Versicolor
    96          2.9           4.2  Versicolor
    97          2.9           4.3  Versicolor
    98          2.5           3.0  Versicolor
    99          2.8           4.1  Versicolor

    [100 rows x 3 columns]

```python
var=[]
for i in range(len(iris_data_12)):
    if iris_data_12['variety'].iloc[i]=='Setosa':
        var.append(1)
    else:
        var.append(-1)
iris_data_12['var'] = var
iris_data_12 = iris_data_12.drop(['variety'],axis=1)
print(iris_data_12)
```

        sepal.width  petal.length  var
    0           3.5           1.4    1
    1           3.0           1.4    1
    2           3.2           1.3    1
    3           3.1           1.5    1
    4           3.6           1.4    1
    ..          ...           ...  ...
    95          3.0           4.2   -1
    96          2.9           4.2   -1
    97          2.9           4.3   -1
    98          2.5           3.0   -1
    99          2.8           4.1   -1

    [100 rows x 3 columns]

    <ipython-input-5-72b04fa0b6fd>:7: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      iris_data_12['var'] = var

```python
iris_data_13 = iris_data[iris_data['variety']!='Versicolor']
print(iris_data_13)
```

         sepal.width  petal.length    variety
    0            3.5           1.4     Setosa
    1            3.0           1.4     Setosa
    2            3.2           1.3     Setosa
    3            3.1           1.5     Setosa
    4            3.6           1.4     Setosa
    ..           ...           ...        ...
    145          3.0           5.2  Virginica
    146          2.5           5.0  Virginica
    147          3.0           5.2  Virginica
    148          3.4           5.4  Virginica
    149          3.0           5.1  Virginica

    [100 rows x 3 columns]

```python
var=[]
for i in range(len(iris_data_13)):
    if iris_data_13['variety'].iloc[i]=='Setosa':
        var.append(1)
    else:
        var.append(-1)
iris_data_13['var'] = var
iris_data_13 = iris_data_13.drop(['variety'],axis=1)
print(iris_data_13)
```

         sepal.width  petal.length  var
    0            3.5           1.4    1
    1            3.0           1.4    1
    2            3.2           1.3    1
    3            3.1           1.5    1
    4            3.6           1.4    1
    ..           ...           ...  ...
    145          3.0           5.2   -1
    146          2.5           5.0   -1
    147          3.0           5.2   -1
    148          3.4           5.4   -1
    149          3.0           5.1   -1

    [100 rows x 3 columns]

    <ipython-input-7-636e89629875>:7: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      iris_data_13['var'] = var

```python
iris_data_23 = iris_data[iris_data['variety']!='Setosa']
print(iris_data_23)
```

         sepal.width  petal.length     variety
    50           3.2           4.7  Versicolor
    51           3.2           4.5  Versicolor
    52           3.1           4.9  Versicolor
    53           2.3           4.0  Versicolor
    54           2.8           4.6  Versicolor
    ..           ...           ...         ...
    145          3.0           5.2   Virginica
    146          2.5           5.0   Virginica
    147          3.0           5.2   Virginica
    148          3.4           5.4   Virginica
    149          3.0           5.1   Virginica

    [100 rows x 3 columns]

```python
var=[]
for i in range(len(iris_data_23)):
    if iris_data_23['variety'].iloc[i]=='Versicolor':
        var.append(1)
    else:
        var.append(-1)
iris_data_23['var'] = var
iris_data_23 = iris_data_23.drop(['variety'],axis=1)
print(iris_data_23)
```

         sepal.width  petal.length  var
    50           3.2           4.7    1
    51           3.2           4.5    1
    52           3.1           4.9    1
    53           2.3           4.0    1
    54           2.8           4.6    1
    ..           ...           ...  ...
    145          3.0           5.2   -1
    146          2.5           5.0   -1
    147          3.0           5.2   -1
    148          3.4           5.4   -1
    149          3.0           5.1   -1

    [100 rows x 3 columns]

    <ipython-input-9-0a34df3c1e97>:7: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame.
    Try using .loc[row_indexer,col_indexer] = value instead

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      iris_data_23['var'] = var

####  Perceptron

```python
data_aug_12 = pd.DataFrame(list(zip(iris_data_12['petal.length'],iris_data_12['sepal.width'], [1 for i in range(len(iris_data_12))],iris_data_12['var'])), columns=['petal.length','sepal.width','augment','var'])
print(data_aug_12)
```

        petal.length  sepal.width  augment  var
    0            1.4          3.5        1    1
    1            1.4          3.0        1    1
    2            1.3          3.2        1    1
    3            1.5          3.1        1    1
    4            1.4          3.6        1    1
    ..           ...          ...      ...  ...
    95           4.2          3.0        1   -1
    96           4.2          2.9        1   -1
    97           4.3          2.9        1   -1
    98           3.0          2.5        1   -1
    99           4.1          2.8        1   -1

    [100 rows x 4 columns]

```python
data_aug_13 = pd.DataFrame(list(zip(iris_data_13['petal.length'],iris_data_13['sepal.width'], [1 for i in range(len(iris_data_13))],iris_data_13['var'])), columns=['petal.length','sepal.width','augment','var'])
print(data_aug_13)
```

        petal.length  sepal.width  augment  var
    0            1.4          3.5        1    1
    1            1.4          3.0        1    1
    2            1.3          3.2        1    1
    3            1.5          3.1        1    1
    4            1.4          3.6        1    1
    ..           ...          ...      ...  ...
    95           5.2          3.0        1   -1
    96           5.0          2.5        1   -1
    97           5.2          3.0        1   -1
    98           5.4          3.4        1   -1
    99           5.1          3.0        1   -1

    [100 rows x 4 columns]

```python
data_aug_23 = pd.DataFrame(list(zip(iris_data_23['petal.length'],iris_data_23['sepal.width'], [1 for i in range(len(iris_data_23))],iris_data_23['var'])), columns=['petal.length','sepal.width','augment','var'])
print(data_aug_23)
```

        petal.length  sepal.width  augment  var
    0            4.7          3.2        1    1
    1            4.5          3.2        1    1
    2            4.9          3.1        1    1
    3            4.0          2.3        1    1
    4            4.6          2.8        1    1
    ..           ...          ...      ...  ...
    95           5.2          3.0        1   -1
    96           5.0          2.5        1   -1
    97           5.2          3.0        1   -1
    98           5.4          3.4        1   -1
    99           5.1          3.0        1   -1

    [100 rows x 4 columns]

**Splitting Features and Labels**

```python
def split(data_aug):
    data_up = data_aug.copy()
    for i in range(len(data_up)):
        if(data_up['var'][i] == 1):
            data_up['petal.length'][i] = - data_up['petal.length'][i]
            data_up['sepal.width'][i] = -data_up['sepal.width'][i]
            data_up['augment'][i] = -data_up['augment'][i]
    # print(data_up)
    data_up = data_up.drop(['var'], axis=1)
    data_up = np.array(data_up)
    # print(data_up)

    labels = data_aug['var']
    labels = np.array(labels)
    # print(labels)

    return (data_up, labels)
```

```python
data_up_12, labels_12 = split(data_aug_12)
data_up_13, labels_13 = split(data_aug_13)
data_up_23, labels_23 = split(data_aug_23)

```

    <ipython-input-13-e04cf10d8a68>:5: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      data_up['petal.length'][i] = - data_up['petal.length'][i]
    <ipython-input-13-e04cf10d8a68>:6: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      data_up['sepal.width'][i] = -data_up['sepal.width'][i]
    <ipython-input-13-e04cf10d8a68>:7: SettingWithCopyWarning:
    A value is trying to be set on a copy of a slice from a DataFrame

    See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
      data_up['augment'][i] = -data_up['augment'][i]

```python
def Perceptron_Train(data_up, labels, eta):

    threshold = 5000
    A = np.zeros(len(data_up[0]))

    # print(A)

    for epoch in range(threshold):
        ctr = 0

        for inp, label in zip(data_up, labels):
            summation = np.dot(inp, A)
            if summation <= 0:
                A += eta*inp
                ctr = 0

            else:
                ctr+=1

        if ctr == len(data_up):     #Until all input is correctly classified
            break
    print(epoch)
    return A

```

```python
x1, x2, b = symbols('x1 x2 1')
P = np.array([x1,x2, b])
print(P)
```

    [x1 x2 1]

```python
eta = 0.01
A_12 = Perceptron_Train(data_up_12, labels_12, eta)
print(A_12)
```

    3
    [ 0.052 -0.041 -0.01 ]

```python
Percep_equation_12 = np.dot(A_12.T, P)
print("Decision Boundary between Setosa and Versicolor: ", Percep_equation_12, ' = 0')
```

    Decision Boundary between Setosa and Versicolor:  -0.01*1 + 0.052*x1 - 0.041*x2  = 0

```python
eta = 0.01
A_13 = Perceptron_Train(data_up_13, labels_13, eta)
print(A_13)
```

    2
    [ 0.032 -0.037 -0.01 ]

```python
Percep_equation_13 = np.dot(A_13.T, P)
print("Decision Boundary between Setosa and Virginica: ", Percep_equation_13, ' = 0')
```

    Decision Boundary between Setosa and Virginica:  -0.01*1 + 0.032*x1 - 0.037*x2  = 0

```python
eta = 0.01
A_23 = Perceptron_Train(data_up_23, labels_23, eta)
print(A_23)
```

    4999
    [ 1.296 -0.837 -3.43 ]

```python
Percep_equation_23 = np.dot(A_23.T, P)
print("Decision Boundary between Versicolor and Virginica: ", Percep_equation_23, ' = 0')
```

    Decision Boundary between Versicolor and Virginica:  -3.42999999999997*1 + 1.29600000000016*x1 - 0.837*x2  = 0

Plotting the Decision Boundary

```python
def plot_percep(A, data_aug):
    a = np.linspace(-1,7,100)
    b = -a*A[0]/A[1] - A[2]/A[1]

    sns.set(style="darkgrid")
    sns.lmplot(x='petal.length',y='sepal.width', data=data_aug, fit_reg=False, hue='var', legend=True)
    plt.plot(a,b)
    plt.show()
```

```python
plot_percep(A_12,data_aug_12)
```

```python
plot_percep(A_23,data_aug_23)
```

```python
plot_percep(A_13, data_aug_13)
```

```python
a = np.linspace(-1,7,100)
b_12 = -a*A_12[0]/A_12[1] - A_12[2]/A_12[1]
b_13 = -a*A_13[0]/A_13[1] - A_13[2]/A_13[1]
b_23 = -a*A_23[0]/A_23[1] - A_23[2]/A_23[1]
sns.set(style="darkgrid")
sns.lmplot(x='petal.length',y='sepal.width', data=iris_data, fit_reg=False, hue='variety', legend=True)
plt.plot(a,b_12)
plt.plot(a,b_13)
plt.plot(a,b_23)
plt.show()
```

####  SVM

```python
def data_svm_split(iris):
    X = iris[['petal.length','sepal.width']]
    X = np.array(X)

    Y = iris['var']
    Y = np.array(Y).reshape(len(iris),1)

    return (X,Y)
```

```python
X_12,Y_12 = data_svm_split(iris_data_12)
X_13,Y_13 = data_svm_split(iris_data_13)
X_23,Y_23 = data_svm_split(iris_data_23)
```

```python
def SVM_Train(X,Y):
    m = len(X)
    P = matrix(np.multiply(np.dot(Y, Y.T), np.dot(X, X.T)))
    q = matrix(np.ones(m) * -1)
    g1 = np.asarray(np.diag(np.ones(m) * -1))

    print(np.array(g1).shape)
    h = matrix(np.zeros(m))
    A = np.reshape((Y.T), (1,m))
    b=[[0]]

    P = matrix(P,(m,m),'d') #dense
    A = matrix(A,(1,m),'d')
    g1 = matrix(g1,(m,m),'d')
    b = matrix(b,(1,1),'d')

    sol = solvers.qp(P, q, g1, h, A, b)
    alpha = np.array(sol['x'])
    ind = (alpha > 1e-4).flatten()
    print(ind)

    W = np.dot(np.transpose(alpha*Y),X)
    print(W)

    for i in range(m):
        if ind[i] == True:
            W0 = Y[i] - np.dot(X[i],W.T)
            print
            break

    print(W0)

    return (W, W0)
```

```python
x1, x2 = symbols('x1 x2')
P_SVM = np.array([x1,x2])
print(P_SVM)
```

    [x1 x2]

```python
W_12,W0_12 = SVM_Train(X_12,Y_12)
```

    (100, 100)
         pcost       dcost       gap    pres   dres
     0: -4.3867e+00 -8.1716e+00  3e+02  1e+01  2e+00
     1: -2.3798e+00 -2.5258e+00  2e+01  1e+00  2e-01
     2: -4.3479e-01 -1.7538e+00  2e+00  5e-02  6e-03
     3: -6.7595e-01 -1.0544e+00  5e-01  1e-02  1e-03
     4: -8.1435e-01 -1.1260e+00  4e-01  4e-03  5e-04
     5: -1.0085e+00 -1.0271e+00  2e-02  8e-05  1e-05
     6: -1.0251e+00 -1.0253e+00  2e-04  9e-07  1e-07
     7: -1.0253e+00 -1.0253e+00  2e-06  9e-09  1e-09
     8: -1.0253e+00 -1.0253e+00  2e-08  9e-11  1e-11
    Optimal solution found.
    [False False False False False False False False False False False False
     False False False False False False False False False False False False
      True False False False False False False False False False False False
     False False False False False  True False False False False False False
     False False False False False False False False False False False False
     False False False False False False False False False False False False
     False False False False False False False False False False False False
     False False False False False False False False False False False False
     False False  True False]
    [[-1.25714286  0.68571428]]
    [1.05714287]

```python
SVM_equation_12 = np.dot(W_12[0].T, P_SVM) + W0_12[0]
print("Decision Boundary between Setosa and Versicolor: ", SVM_equation_12, ' = 0')
```

    Decision Boundary between Setosa and Versicolor:  -1.25714285940527*x1 + 0.685714284290053*x2 + 1.05714286628384  = 0

```python
W_13,W0_13 = SVM_Train(X_13,Y_13)
```

    (100, 100)
         pcost       dcost       gap    pres   dres
     0: -3.3137e+00 -4.9765e+00  2e+02  2e+01  1e+00
     1: -1.6126e-01 -8.1531e-01  1e+01  9e-01  8e-02
     2:  1.9106e-02 -6.5205e-01  9e-01  1e-02  1e-03
     3: -1.5487e-01 -2.7974e-01  1e-01  7e-04  6e-05
     4: -2.2572e-01 -2.8972e-01  6e-02  1e-04  1e-05
     5: -2.5983e-01 -2.6459e-01  5e-03  1e-05  8e-07
     6: -2.6415e-01 -2.6421e-01  5e-05  1e-07  9e-09
     7: -2.6420e-01 -2.6420e-01  5e-07  1e-09  9e-11
     8: -2.6420e-01 -2.6420e-01  5e-09  1e-11  9e-13
    Optimal solution found.
    [False False False False False False False False False False False False
     False False False False False False False False False False False False
      True False False False False False False False False False False False
     False False False False False False False False False False False False
     False False False False False False False False  True False False False
     False False False False False False False False False False False False
     False False False False False False False False False False False False
     False False False False False False False False False False False False
     False False False False]
    [[-0.68692206  0.23778071]]
    [1.49669752]

```python
SVM_equation_13 = np.dot(W_13[0].T, P_SVM) + W0_13[0]
print("Decision Boundary between Setosa and Virginica: ", SVM_equation_13, ' = 0')
```

    Decision Boundary between Setosa and Virginica:  -0.686922063593103*x1 + 0.237780707125038*x2 + 1.49669751660177  = 0

```python
W_23,W0_23 = SVM_Train(X_23,Y_23)
```

    (100, 100)
         pcost       dcost       gap    pres   dres
     0: -3.7394e+01 -8.7261e+01  5e+02  2e+01  2e+00
     1: -1.0739e+02 -1.7754e+02  3e+02  1e+01  2e+00
     2: -5.7593e+02 -8.7585e+02  4e+02  1e+01  2e+00
     3: -1.9989e+03 -2.2838e+03  3e+02  1e+01  1e+00
     4: -6.2946e+03 -6.9640e+03  7e+02  1e+01  1e+00
     5: -3.5722e+04 -3.7828e+04  2e+03  1e+01  1e+00
     6: -1.6095e+05 -1.6859e+05  8e+03  1e+01  1e+00
     7: -2.7868e+06 -2.8212e+06  3e+04  9e+00  1e+00
     8: -1.2159e+08 -1.2231e+08  7e+05  9e+00  1e+00
     9: -1.9001e+10 -1.9011e+10  1e+07  9e+00  1e+00
    10: -1.9282e+10 -1.9292e+10  1e+07  9e+00  1e+00
    11: -1.9287e+10 -1.9297e+10  1e+07  9e+00  1e+00
    12: -2.2196e+10 -2.2207e+10  1e+07  9e+00  1e+00
    13: -3.9214e+10 -3.9235e+10  2e+07  9e+00  1e+00
    14: -7.6312e+10 -7.6348e+10  4e+07  9e+00  1e+00
    15: -9.6129e+10 -9.6168e+10  4e+07  9e+00  1e+00
    Terminated (singular KKT matrix).
    [ True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True  True  True  True  True  True  True  True  True
      True  True  True  True]
    [[ 0.00000000e+00 -3.05175781e-05]]
    [1.00009766]

```python
SVM_equation_23 = np.dot(W_23[0].T, P_SVM) + W0_23[0]
print("Decision Boundary between Versicolor and Virginica: ", SVM_equation_23, ' = 0')
```

    Decision Boundary between Versicolor and Virginica:  1.00009765625 - 3.0517578125e-5*x2  = 0

```python
def plot_SVM(W,W0, data_aug):
    a = np.linspace(-1,7,100)
    b = -a*W[0][0]/W[0][1] - W0[0]/W[0][1]
    sns.set(style="darkgrid")
    sns.lmplot(x='petal.length',y='sepal.width', data=data_aug, fit_reg=False, hue='var', legend=True)
    plt.plot(a,b)
    plt.show()
```

```python
plot_SVM(W_12,W0_12,iris_data_12)
```

```python
plot_SVM(W_13,W0_13,iris_data_13)
```

```python
plot_SVM(W_23,W0_23,iris_data_23)
```

```python
a = np.linspace(-1,7,100)
b_12 = -a*W_12[0][0]/W_12[0][1] - W0_12[0]/W_12[0][1]
b_13 = -a*W_13[0][0]/W_13[0][1] - W0_13[0]/W_13[0][1]
sns.set(style="darkgrid")
sns.lmplot(x='petal.length',y='sepal.width', data=iris_data, fit_reg=False, hue='variety', legend=True)

plt.plot(a,b_12)
plt.plot(a,b_13)
plt.show()
```

```python
a = np.linspace(-1,7,100)
b_12 = -a*W_12[0][0]/W_12[0][1] - W0_12[0]/W_12[0][1]
b_13 = -a*W_13[0][0]/W_13[0][1] - W0_13[0]/W_13[0][1]
b_23 = -a*W_23[0][0]/W_23[0][1] - W0_23[0]/W_23[0][1]
sns.set(style="darkgrid")
sns.lmplot(x='petal.length',y='sepal.width', data=iris_data, fit_reg=False, hue='variety', legend=True)
plt.plot(a,b_12)
plt.plot(a,b_13)
plt.plot(a,b_23)
plt.show()
```

max_stars_repo_name: PaletiKrishnasai/Pattern-Recognition

max_issues_repo_name: PaletiKrishnasai/Pattern-Recognition

max_forks_repo_name: PaletiKrishnasai/Pattern-Recognition