Note

Go to the end to download the full example code.

Dustywave TVA dispersion relation#

This example shows how to derive the dustywave TVA dispersion relation using SymPy

 8 import matplotlib.pyplot as plt
 9 import numpy as np
10 import sympy as sp

Usefull symbols

14 omega = sp.symbols(r"\omega", complex=True)
15
16 k, cs, ts, eps = sp.symbols(
17     r"k c_s t_s \epsilon",
18     positive=True,
19     real=True,
20 )
21
22 i = sp.I  # imaginary unit
23
24 a = k * cs  # s^-1
25 b = a**2 * ts * eps  # s^-1

The perturbation matrix is M = K(omega = 0) The eigenfrequencies are the roots of det(M + i omega I) = 0

31 K = sp.Matrix(
32     [
33         [i * omega, 0, -i * a],
34         [b * (1 - eps), i * omega - b, 0],
35         [-i * a * (1 - eps), i * a, i * omega],
36     ]
37 )  # s^-1

Compute the determinant

42 det = sp.factor(K.det())
43 print(sp.latex(det))

\omega \left(\epsilon \omega c_{s}^{2} k^{2} t_{s} - i \epsilon c_{s}^{2} k^{2} - i \omega^{2} + i c_{s}^{2} k^{2}\right)

Let’s remove the leading omega mode (omega = 0)

47 det /= i * omega
48 det = sp.simplify(det)
49 det = sp.collect(det, omega)
50 print(sp.latex(det))

- i \epsilon \omega c_{s}^{2} k^{2} t_{s} - \epsilon c_{s}^{2} k^{2} - \omega^{2} + c_{s}^{2} k^{2}

Find the roots of the dispersion relation

54 r1, r2 = sp.solve(sp.Eq(det, 0), omega)
55 print(sp.latex(r1))
56 print(sp.latex(r2))
57
58 print(r1)
59 print(r2)

\frac{c_{s} k \left(- i \epsilon c_{s} k t_{s} - \sqrt{- \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 4 \epsilon + 4}\right)}{2}
\frac{c_{s} k \left(- i \epsilon c_{s} k t_{s} + \sqrt{- \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 4 \epsilon + 4}\right)}{2}
c_s*k*(-I*\epsilon*c_s*k*t_s - sqrt(-\epsilon**2*c_s**2*k**2*t_s**2 - 4*\epsilon + 4))/2
c_s*k*(-I*\epsilon*c_s*k*t_s + sqrt(-\epsilon**2*c_s**2*k**2*t_s**2 - 4*\epsilon + 4))/2

Function to plot the roots

 64 def get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value):
 65
 66     # Substitute all parameters
 67     r1_num = _r1.subs({cs: cs_value, ts: ts_value, eps: eps_value})
 68
 69     r2_num = _r2.subs({cs: cs_value, ts: ts_value, eps: eps_value})
 70
 71     r1_num_re = sp.re(r1_num)
 72     r1_num_im = sp.im(r1_num)
 73     r2_num_re = sp.re(r2_num)
 74     r2_num_im = sp.im(r2_num)
 75
 76     # Lambdify only k remains
 77     r1_re_func = sp.lambdify(k, r1_num_re, modules="numpy")
 78     r1_im_func = sp.lambdify(k, r1_num_im, modules="numpy")
 79     r2_re_func = sp.lambdify(k, r2_num_re, modules="numpy")
 80     r2_im_func = sp.lambdify(k, r2_num_im, modules="numpy")
 81
 82     # Evaluate
 83     r1_vals_re = r1_re_func(k_list)
 84     r1_vals_im = r1_im_func(k_list)
 85     r2_vals_re = r2_re_func(k_list)
 86     r2_vals_im = r2_im_func(k_list)
 87
 88     def restore(lst):
 89         # if it is not a numpy array return a np.zeros_like(k_list)
 90         if not isinstance(lst, np.ndarray):
 91             return np.zeros_like(k_list)
 92         return lst
 93
 94     r1_vals_re = restore(r1_vals_re)
 95     r1_vals_im = restore(r1_vals_im)
 96     r2_vals_re = restore(r2_vals_re)
 97     r2_vals_im = restore(r2_vals_im)
 98
 99     return r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im
100
101
102 def get_roots_LP14(k_list, eps_value, cs_value, ts_value):
103     _r1 = +cs * sp.sqrt(1 - eps) * k - i * ts * k**2 * cs**2 * eps / 2
104     _r2 = -cs * sp.sqrt(1 - eps) * k - i * ts * k**2 * cs**2 * eps / 2
105
106     return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value)
107
108
109 def get_overroots_DCL26_simple(k_list, eps_value, cs_value, ts_value):
110     _r1 = +cs * sp.sqrt(1 - eps) * k + i * k**2 * cs**2 * eps * ts * (-1 + 1) / 2
111     _r2 = -cs * sp.sqrt(1 - eps) * k + i * k**2 * cs**2 * eps * ts * (-1 - 1) / 2
112
113     return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value)
114
115
116 def get_overroots_DCL26(k_list, eps_value, cs_value, ts_value):
117     D = 4 * (1 - eps) - eps**2 * cs**2 * ts**2 * k**2
118
119     sqrtD_real = sp.sqrt(sp.Max(D, 0))
120     sqrtD_imag = sp.sqrt(sp.Max(-D, 0))
121
122     _r1 = cs * k / 2 * (+sqrtD_real + i * (sqrtD_imag - eps * cs * k * ts))
123
124     _r2 = cs * k / 2 * (-sqrtD_real + i * (-sqrtD_imag - eps * cs * k * ts))
125
126     print(sp.latex(sp.Abs(_r1)))
127     print(sp.latex(sp.Abs(_r2)))
128
129     return get_roots(_r1, _r2, k_list, eps_value, cs_value, ts_value)
130
131
132 def plot_case(k_plot, eps_value, cs_value, ts_value):
133
134     r1_vals_re, r1_vals_im, r2_vals_re, r2_vals_im = get_roots(
135         r1, r2, k_plot, eps_value, cs_value, ts_value
136     )
137
138     r1_vals_re_LP14, r1_vals_im_LP14, r2_vals_re_LP14, r2_vals_im_LP14 = get_roots_LP14(
139         k_plot, eps_value, cs_value, ts_value
140     )
141
142     (
143         r1_vals_re_DCL26_simple,
144         r1_vals_im_DCL26_simple,
145         r2_vals_re_DCL26_simple,
146         r2_vals_im_DCL26_simple,
147     ) = get_overroots_DCL26_simple(k_plot, eps_value, cs_value, ts_value)
148     r1_vals_re_DCL26, r1_vals_im_DCL26, r2_vals_re_DCL26, r2_vals_im_DCL26 = get_overroots_DCL26(
149         k_plot, eps_value, cs_value, ts_value
150     )
151
152     # Create figure
153     fig, axs = plt.subplots(2, 2, figsize=(8, 8), sharex=True)
154
155     # Real parts
156     axs[0, 0].plot(k_plot, r1_vals_re, color="0", linewidth=2, label=r"Re($\omega_+$)")
157     axs[0, 0].plot(k_plot, r2_vals_re, color="0", linewidth=2, label=r"Re($\omega_-$)")
158     axs[0, 0].plot(k_plot, r1_vals_re_LP14, "--", label=r"Re($\omega_{+,LP14}$)")
159     axs[0, 0].plot(k_plot, r2_vals_re_LP14, "--", label=r"Re($\omega_{-,LP14}$)")
160     axs[0, 0].plot(
161         k_plot, r1_vals_re_DCL26_simple, linestyle="dotted", label=r"Re($\omega_{+,approx}$)"
162     )
163     axs[0, 0].plot(
164         k_plot, r2_vals_re_DCL26_simple, linestyle="dotted", label=r"Re($\omega_{-,approx}$)"
165     )
166     # axs[0,0].plot(k_plot, r1_vals_re_DCL26,"--", label="Re($r_1$) DCL26")
167     # axs[0,0].plot(k_plot, r2_vals_re_DCL26,"--", label="Re($r_2$) DCL26")
168     axs[0, 0].set_ylabel("Real part")
169     axs[0, 0].grid(True)
170     axs[0, 0].legend()
171
172     # Imaginary parts
173     axs[0, 1].plot(k_plot, r1_vals_im, color="0", linewidth=2, label=r"Im($\omega_+$)")
174     axs[0, 1].plot(k_plot, r2_vals_im, color="0", linewidth=2, label=r"Im($\omega_-$)")
175     axs[0, 1].plot(k_plot, r1_vals_im_LP14, "--", label=r"Im($\omega_{+,LP14}$)")
176     axs[0, 1].plot(k_plot, r2_vals_im_LP14, "--", label=r"Im($\omega_{-,LP14}$)")
177     axs[0, 1].plot(
178         k_plot, r1_vals_im_DCL26_simple, linestyle="dotted", label=r"Im($\omega_{+,approx}$)"
179     )
180     axs[0, 1].plot(
181         k_plot, r2_vals_im_DCL26_simple, linestyle="dotted", label=r"Im($\omega_{-,approx}$)"
182     )
183     # axs[0,1].plot(k_plot, r1_vals_im_DCL26,"--", label="Im($r_1$) DCL26")
184     # axs[0,1].plot(k_plot, r2_vals_im_DCL26,"--", label="Im($r_2$) DCL26")
185     axs[0, 1].set_xlabel("$k$")
186     axs[0, 1].set_ylabel("Imaginary part")
187     axs[0, 1].grid(True)
188     axs[0, 1].legend()
189
190     # Abs
191     r1_vals_abs = np.sqrt(r1_vals_re**2 + r1_vals_im**2)
192     r2_vals_abs = np.sqrt(r2_vals_re**2 + r2_vals_im**2)
193     r1_vals_abs_LP14 = np.sqrt(r1_vals_re_LP14**2 + r1_vals_im_LP14**2)
194     r2_vals_abs_LP14 = np.sqrt(r2_vals_re_LP14**2 + r2_vals_im_LP14**2)
195     r1_vals_abs_DCL26_simple = np.sqrt(r1_vals_re_DCL26_simple**2 + r1_vals_im_DCL26_simple**2)
196     r2_vals_abs_DCL26_simple = np.sqrt(r2_vals_re_DCL26_simple**2 + r2_vals_im_DCL26_simple**2)
197     r1_vals_abs_DCL26 = np.sqrt(r1_vals_re_DCL26**2 + r1_vals_im_DCL26**2)
198     r2_vals_abs_DCL26 = np.sqrt(r2_vals_re_DCL26**2 + r2_vals_im_DCL26**2)
199     axs[1, 0].plot(k_plot, r1_vals_abs, color="0", linewidth=2, label=r"Abs($\omega_+$)")
200     axs[1, 0].plot(k_plot, r2_vals_abs, color="0", linewidth=2, label=r"Abs($\omega_-$)")
201     axs[1, 0].plot(k_plot, r1_vals_abs_LP14, "--", label=r"Abs($\omega_{+,LP14}$)")
202     axs[1, 0].plot(k_plot, r2_vals_abs_LP14, "--", label=r"Abs($\omega_{-,LP14}$)")
203     axs[1, 0].plot(
204         k_plot, r1_vals_abs_DCL26_simple, linestyle="dotted", label=r"Abs($\omega_{+,approx}$)"
205     )
206     axs[1, 0].plot(
207         k_plot, r2_vals_abs_DCL26_simple, linestyle="dotted", label=r"Abs($\omega_{-,approx}$)"
208     )
209
210     def approx(_k, _cs, _ts, _eps):
211         print(type(_k), type(_cs), type(_ts), type(_eps))
212         return _cs * _k * np.sqrt((1 - _eps) + (_k * _cs * _ts * _eps) ** 2)
213
214     axs[1, 0].plot(
215         k_plot,
216         approx(k_plot, cs_value, ts_value, eps_value),
217         "--",
218         label=r"$max(\vert \omega_{\pm,approx} \vert)$",
219     )
220
221     # axs[1,0].plot(k_plot, r1_vals_abs_DCL26,"--", label="Abs($r_1$) DCL26")
222     # axs[1,0].plot(k_plot, r2_vals_abs_DCL26,"--", label="Abs($r_2$) DCL26")
223     axs[1, 0].set_xlabel("$k$")
224     axs[1, 0].set_ylabel("Abs part")
225     axs[1, 0].grid(True)
226     axs[1, 0].legend()
227
228     # delta with max
229     r_max = np.maximum(r1_vals_abs, r2_vals_abs)
230     r_max_LP14 = np.maximum(r1_vals_abs_LP14, r2_vals_abs_LP14)
231     r_max_DCL26_simple = np.maximum(r1_vals_abs_DCL26_simple, r2_vals_abs_DCL26_simple)
232     r_max_DCL26 = np.maximum(r1_vals_abs_DCL26, r2_vals_abs_DCL26)
233     axs[1, 1].plot(k_plot, (r_max_LP14 - r_max) / r_max, label="Ana - LP14")
234     axs[1, 1].plot(k_plot, (r_max_DCL26_simple - r_max) / r_max, label="Ana - DCL26 simple")
235     # axs[1,1].plot(k_plot, (r_max_DCL26 - r_max) / r_max, label="Ana - DCL26")
236     axs[1, 1].set_xlabel("$k$")
237     axs[1, 1].set_ylabel("Abs(Ana) - Abs(Max model) / Abs(Ana)")
238     axs[1, 1].grid(True)
239     axs[1, 1].legend()
240
241     plt.suptitle(f"eps = {eps_value}, cs = {cs_value}, ts = {ts_value}")
242
243     plt.tight_layout()

Plot the case eps = 0.5, cs = 1.0, ts = 1.0

248 k_plot = np.linspace(0, 5, 1000)
249 plot_case(k_plot, 0.5, 1.0, 1.0)
250 plt.show()
eps = 0.5, cs = 1.0, ts = 1.0
\frac{c_{s} k \sqrt{\epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 2 \epsilon c_{s} k t_{s} \sqrt{\max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)} + \max\left(0, - \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 4 \epsilon + 4\right) + \max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)}}{2}
\frac{c_{s} k \sqrt{\epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 2 \epsilon c_{s} k t_{s} \sqrt{\max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)} + \max\left(0, - \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 4 \epsilon + 4\right) + \max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)}}{2}
<class 'numpy.ndarray'> <class 'float'> <class 'float'> <class 'float'>
/work/doc/sphinx/examples/physics/run_dustywave_sympy.py:233: RuntimeWarning: invalid value encountered in divide
  axs[1, 1].plot(k_plot, (r_max_LP14 - r_max) / r_max, label="Ana - LP14")
/work/doc/sphinx/examples/physics/run_dustywave_sympy.py:234: RuntimeWarning: invalid value encountered in divide
  axs[1, 1].plot(k_plot, (r_max_DCL26_simple - r_max) / r_max, label="Ana - DCL26 simple")

Plot the case eps = 0.1, cs = 1.0, ts = 1.0

255 k_plot = np.linspace(0, 40, 1000)
256 plot_case(k_plot, 0.1, 1.0, 1.0)
257 plt.show()
eps = 0.1, cs = 1.0, ts = 1.0
\frac{c_{s} k \sqrt{\epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 2 \epsilon c_{s} k t_{s} \sqrt{\max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)} + \max\left(0, - \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 4 \epsilon + 4\right) + \max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)}}{2}
\frac{c_{s} k \sqrt{\epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 2 \epsilon c_{s} k t_{s} \sqrt{\max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)} + \max\left(0, - \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} - 4 \epsilon + 4\right) + \max\left(0, \epsilon^{2} c_{s}^{2} k^{2} t_{s}^{2} + 4 \epsilon - 4\right)}}{2}
<class 'numpy.ndarray'> <class 'float'> <class 'float'> <class 'float'>

Total running time of the script: (0 minutes 3.917 seconds)

Estimated memory usage: 160 MB

Gallery generated by Sphinx-Gallery