```python
from sympy import *
init_printing(use_latex='mathjax')
Re,Da,theta_r,D,rho,rho_ph,L_x,lam,tau,dt,k,x = symbols('Re Da theta_r D rho rho_ph L_x lambda tau dt 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
dt  = solve(Da - L_x**2*rho_ph/D/dt,dt)[0]
```


```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) + rho_ph*H/dt # 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 Da - \frac{Da \lambda}{Re} - Re \lambda + \pi^{2} k^{2} + \lambda^{2}$




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




$\displaystyle \sqrt{Da + \pi^{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,Da],lam,"numpy")
Ref       = lambdify(Da,Re_opt,"numpy")
Re1       = np.linspace(np.pi/2,5*np.pi,500)
Da1       = np.linspace(0,100,500)
[Re2,Da2] = np.meshgrid(Re1,Da1)
lam2      = np.stack(lamf(Re2,Da2),axis=0)
lam2      = np.nanmin(lam2,axis=0)
Re_opt1   = Ref(Da1)
lam1      = np.stack(lamf(Re_opt1,Da1),axis=0)
lam1      = np.nanmin(lam1,axis=0)
```


```python
%matplotlib notebook
from matplotlib import pyplot as plt
from matplotlib import cm
plt.rcParams.update({"text.usetex": True, "font.size": 14})
f   = plt.figure(figsize=(7,5))
ax  = plt.subplot()
pos = ax.contourf(Re2,Da2,lam2,10,cmap=cm.cool)
ax.contour(Re2,Da2,lam2,10,linewidths=1,colors='k')
ax.plot(Re_opt1,Da1,ls='--',linewidth=4,color='w')
ax.set_xlabel("$Re$"); ax.set_ylabel("$Da$")
ax.set_xticks([np.pi/2,np.pi,2*np.pi,3*np.pi,4*np.pi,5*np.pi])
ax.set_xticklabels(["$\pi/2$","$\pi$","$2\pi$","$3\pi$","4$\pi$","5$\pi$"])
ax.set_aspect(1.0/ax.get_data_ratio(), adjustable='box')
cbar = f.colorbar(pos,ax=ax)
cbar.set_label("$\mathrm{min}\{\Re(\lambda_k)\}$")
```


    <IPython.core.display.Javascript object>







```python

```

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