Loading [MathJax]/jax/output/CommonHTML/fonts/TeX/fontdata.js

%matplotlib inline

Rownanie adwekcji

Dla przepływu jednowymiarowego ze stałą prędkością v = c, równanie to przyjmuje postać:

$$\frac{\partial \rho}{\partial t} + c \frac{\partial \rho}{\partial x} = 0$$

Równanie to opisuje unoszenie pewnej substancji przez przepływający płyn. Nietrudno sprawdzić, że dowolna funkcja w postaci fali biegnącej: $u(x, t) = f (x − ct)$, spełnia to równanie.

Metoda Laxa-Wendroffa

$$u_{i,j+1}=\frac{1}{2}\beta(1+\beta)u_{i-1,j}+(1-\beta^2)u_{i,j}-\frac{1}{2}\beta(1-\beta)u_{i+1,j}$$

import numpy as np; import matplotlib.pyplot as plt;
from mpl_toolkits import mplot3d
from matplotlib import animation
from IPython.display import HTML

m = 100; c = 1.
dx = 1./m; x = np.linspace(0,1,m+1); 
beta = 0.8
# beta = c*dt / dx
dt = beta*dx / c 
Tfinal = 0.5
n = int(Tfinal/dt); t = np.linspace(0,Tfinal,n+1)
u = np.zeros((n+1, m+1), float); u[0] = np.exp(-300.*(x-0.2)**2)

for j in range(0, n):
    for i in range(0, m-1):
        u[j+1][i+1] = (1. - beta*beta)*u[j][i+1]-(0.5*beta)*(1.-beta)*u[j][i+2] +(0.5*beta)*(1. + beta)*u[j][i]

X, T = np.meshgrid(x, t)
fig=plt.figure(figsize=(10,8), dpi= 100, facecolor='w', edgecolor='k')
ax = plt.axes(projection ='3d')
ax.plot_wireframe(X, T, u, cmap ='viridis') 
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')

/home/michal/.local/lib/python3.6/site-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)

Text(0.5, 0, 'u')

%%capture

fig, ax = plt.subplots()


line, = ax.plot(x, u[0])


def animate(i):
    line.set_ydata(u[i])  # update the data.
    return line,


ani = animation.FuncAnimation(
    fig, animate, interval=1000/n, blit=True, save_count=n)

HTML(ani.to_jshtml())

Once Loop Reflect

Rownanie Burgersa

$$\frac{\partial u}{\partial t} + \epsilon \frac{\partial (u^2/2)}{\partial x} = 0$$

Mozna na nie spojrzeć jak na równanie adwekcji, w którym prędkość fali c = $\epsilon$u jest proporcjonalna do amplitudy. Powoduje to zmianę kształtu fali w czasie. Wyższe części fali będą przemiaeszczać się na front, w efekcie tworząc ostry brdzeg, nazywany falą uderzeniową.

Algorytm skokowy

$$u_{i,j+1} = u_{i,j-1} - \beta \left[\frac{u_{i+1,j}^2-u_{i-1,j}^2}{2}\right], \beta=\frac{\epsilon}{\Delta x / \Delta t}$$

m = 100; dx = 2*np.pi/m; x = np.linspace(0,2*np.pi,m+1); 
beta = 0.2; eps = 1; dt = beta*dx/eps; T_final = 0.8;
n = int(T_final/dt); t = np.linspace(0,T_final,n+1)

u0 = 2*np.sin(x); u=np.zeros((n+1, m+1), float)

u[0]=u0; u[1]=u0

for j in range(1,n):
    for i in range(1,m):
        u[j+1][i]=u[j-1][i]-beta*0.5*((u[j][i+1])**2-(u[j][i-1])**2)

X, T = np.meshgrid(x, t)

fig=plt.figure(figsize=(10,8), dpi= 100, facecolor='w', edgecolor='k')
ax = plt.axes(projection ='3d')
ax.plot_wireframe(X, T, u, cmap ='viridis') 
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')

/home/michal/.local/lib/python3.6/site-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)

Text(0.5, 0, 'u')

%%capture

fig, ax = plt.subplots()


line, = ax.plot(x, u[0])


def animate(i):
    line.set_ydata(u[i])  # update the data.
    return line,


ani = animation.FuncAnimation(
    fig, animate, interval=1000/n, blit=True, save_count=n)

HTML(ani.to_jshtml())

Once Loop Reflect

Metoda Laxa-Wedroffa

$$u_{i,j+1} = u_{i,j} - \frac{\beta}{4} \left(u_{i+1,j}^2-u_{i-1,j}^2\right) +\frac{\beta^2}{8}\left[(u_{i+1,j}+u_{i,j})\left(u_{i+1,j}^2-u_{i,j}^2\right)-(u_{i,j}+u_{i-1,j})\left(u_{i,j}^2-u_{i-1,j}^2\right)\right] , \beta=\frac{\epsilon}{\Delta x / \Delta t}$$

w0 = 2*np.sin(x); w=np.zeros((n+1, m+1), float)

w[0]=w0

for j in range(0,n):
    for i in range(1,m):
        w[j+1][i]=w[j][i]-beta*0.25*((w[j][i+1])**2-(w[j][i-1])**2)
        +(beta**2)*0.125*((w[j][i+1]+w[j][i])*((w[j][i+1])**2-(w[j][i])**2)
                         -(w[j][i]+w[j][i-1])*((w[j][i])**2-(w[j][i-1])**2))

fig=plt.figure(figsize=(10,8), dpi= 100, facecolor='w', edgecolor='k')
ax = plt.axes(projection ='3d')
ax.plot_wireframe(X, T, w, cmap ='viridis') 
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_zlabel('u')

/home/michal/.local/lib/python3.6/site-packages/numpy/core/_asarray.py:136: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray
  return array(a, dtype, copy=False, order=order, subok=True)

Text(0.5, 0, 'u')

%%capture

fig, ax = plt.subplots()


line, = ax.plot(x, w[0])


def animate(i):
    line.set_ydata(w[i])  # update the data.
    return line,


ani = animation.FuncAnimation(
    fig, animate, interval=1000/n, blit=True, save_count=n);

HTML(ani.to_jshtml())

Once Loop Reflect