```python
from sympy import *
init_printing(use_latex='mathjax')
Re,theta_r,D,rho,L_x,lam,tau,k,x = symbols('Re theta_r D rho L_x lambda tau k x', positive=True)
C0 = symbols('C0')
```


```python
rho = solve(Re - rho*sqrt(D/rho/theta_r)*L_x/D,rho)[0]  # density from Reynolds number Re
V_p = sqrt(D/rho/theta_r)                               # velocity scale - wave velocity
```


```python
H      = C0*exp(-lam*tau*V_p/L_x)*sin(pi*k*x/L_x)
q      = C0*exp(-lam*tau*V_p/L_x)*cos(pi*k*x/L_x)
```


```python
eq1 =     rho*diff(H,tau)     +   diff(q,x)             # mass balance
eq2 = theta_r*diff(q,tau) + q + D*diff(H,x)             # momentum balance
```


```python
disp = diff(eq1,tau)*theta_r - diff(eq2,x) + eq1
disp = expand(disp)
disp = expand(simplify(disp/disp.coeff(lam**2)))
disp
```




$\displaystyle - Re \lambda + \pi^{2} k^{2} + \lambda^{2}$




```python
sol     = solve(disp,lam)
disc    = discriminant(disp,lam)
Re_opt  = solve(disc,Re)[0].subs(k,1)
Re_opt
```




$\displaystyle 2 \pi$




```python
from sympy import maximum, lambdify
import numpy as np
import warnings
warnings.filterwarnings('ignore')
lam  = [re(s.subs(k,1)) for s in sol]
lamf = lambdify(Re,lam,"numpy")
Re1  = np.linspace(0,2*float(Re_opt),1000)
lam1 = np.stack(lamf(Re1),axis=0)
lam1 = np.nanmin(lam1,axis=0)
```


```python
%matplotlib notebook
from matplotlib import pyplot as plt
plt.rcParams.update({"text.usetex": True, "font.size": 14})
f  = plt.figure(figsize=(5,5))
ax = plt.subplot()
ax.plot(Re1,lam1)
ax.axvline(x=Re_opt.evalf(),ls='--',c='k',linewidth=2)
ax.set_xlabel("$Re$"); ax.set_ylabel("$\mathrm{min}\{\Re(\lambda_k)\}$")
ax.set_xticks([0,np.pi,2*np.pi,3*np.pi,4*np.pi])
ax.set_xticklabels(["$0$","$\pi$","$2\pi$","$3\pi$","$4\pi$"])
ax.set_xlim(0,4*np.pi)
ax.grid()
```


    <IPython.core.display.Javascript object>







```python

```

c6c9f411c1bf47625bd7e1f8312eeed8708ba627
ipynb
Jupyter Notebook
dispersion_analysis/dispersion_analysis_stationary_diffusion1D.ipynb
PTsolvers/PseudoTransientDiffusion.jl
07b3e2e52d04a3f3f9e8fb724bca740ef57249df
MIT
2021-12-06T19:25:10.000Z
2021-12-06T19:25:10.000Z
dispersion_analysis/dispersion_analysis_stationary_diffusion1D.ipynb
PTsolvers/PseudoTransientDiffusion.jl
07b3e2e52d04a3f3f9e8fb724bca740ef57249df
MIT
dispersion_analysis/dispersion_analysis_stationary_diffusion1D.ipynb
PTsolvers/PseudoTransientDiffusion.jl
07b3e2e52d04a3f3f9e8fb724bca740ef57249df
MIT
Qwen/Qwen-72B
1. YES
2. YES