Численное решение уравнения волны Зоммерфельда в проводе

[palz]

Здравствуйте!
Пытаюсь численно решить уравнение для того, чтобы получить постоянную распространения волны Зоммерфельда $\beta$ в проводе.
Уравнение имеет вид $1+\frac{\varepsilon_mk}{k_m}\frac{I_1(k_ma)K_0(ka)}{I_0(k_ma)K_1(ka)}=0$. Здесь $ k=(\beta^2-k_0^2)^{1/2}$, $k_0=\omega/c $, $ k_m=(\beta^2-\varepsilon_mk_0^2)^{1/2}$, $\varepsilon_m=1-\frac{\omega_p^2}{\omega(\omega+i\gamma)}$ - диэлектрическая проницаемость материала провода, a - радиус провода, I и K - модифицированные функции Бесселя.
Значения параметров
$\omega_p=2\cdot10^{13} c^{-1}$
$\gamma=2,5\cdot10^{12} c^{-1}$
a=100 мкм
Необходимо найти $\beta$ в интервале частот $\nu$ от 0,5 до 4 ТГц.
Пытаюсь решить численно в пакете Mathematica встроенной функцией FindRoot, задав начальное приближение для волны вдоль плоской поверхности (пов. плазмон) $\beta=k_0(\varepsilon_m/(\varepsilon_m+1))^{1/2}$, однако на частотах больших 2.9 ТГц вываливается "FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances."
Пытался также предварительно провести минимизацию функцией NMinimize, пробовал определять функции Бесселя и Макдональда как функции, вычисляемые с точностью до 100 знаков, варьировать границы поиска решения, умножать уравнение на знаменатель, увеличивать количество итераций в функциях NMinimize и FindRoot. Все это не дает результатов.

[Pphantom]

Выберите систему единиц для задачи так, чтобы все (или хотя бы большинство) расчетных величин были порядка единицы.

[Vince Diesel]

Приведите весь код, будет больше шансов получить содержательный ответ.

[amon]

palz в сообщении #1421405 писал(а):

Уравнение имеет вид $1+\frac{\varepsilon_mk}{k_m}\frac{I_1(k_ma)K_0(ka)}{I_0(k_ma)K_1(ka)}=0$.

Вы в этом уверены? Насколько я помню, должно быть что-то вроде $1-\frac{\varepsilon_mk}{k_m}\frac{I_1(k_ma)K_0(ka)}{I_0(k_ma)K_1(ka)}=0$.

[palz]

Перешел к новым единицам. Теперь частота измеряется в ТГц, расстояния - в сотнях мкм, соответственно скорость света - c = 2.99792458.
Вот код

\[Nu]min = 0.5; \[Nu]max = 4.0;
N\[Nu] = 20;
\[Nu]a = Array[# &, N\[Nu], {\[Nu]min, \[Nu]max}]; \[Omega]a = 
 2*\[Pi]*\[Nu]a;
\[Beta]a = Array[0 &, N\[Nu]]; k0a = \[Omega]a/c;
c = 2.99792458;
a = 1.0;
\[Omega]p = 2*10;
\[Gamma] = 2.5*\[Omega]p/20;
\[Epsilon]m[\[Omega]_] := 
  1 - \[Omega]p^2/(\[Omega]*(\[Omega] + I*\[Gamma]));
f[\[Omega]_] := (\[Omega]/c)*
   Sqrt[\[Epsilon]m[\[Omega]]/(\[Epsilon]m[\[Omega]] + 1)];
f1[\[Beta]r_, \[Beta]i_, \[Omega]_] := 
  Block[{kp = Sqrt[(\[Beta]r + I*\[Beta]i)^2 - (\[Omega] /c)^2], 
    km = Sqrt[(\[Beta]r + 
         I*\[Beta]i)^2 - \[Epsilon]m[\[Omega]]*(\[Omega] /
          c)^2]}, (1 + \[Epsilon]m[\[Omega]]*
      Sqrt[((\[Beta]r + I*\[Beta]i)^2 - (\[Omega]/c)^2)/((\[Beta]r + 
          I*\[Beta]i)^2 - \[Epsilon]m[\[Omega]]*(\[Omega]/c)^2)]*
      BesselI[1, km*a]/BesselI[0, km*a]*BesselK[0, kp*a]/
      BesselK[1, kp*a])
   ];
fa = f /@ \[Omega]a;
For[s = 1, s <= N\[Nu], s++,
  {\[Beta]r1, \[Beta]i1} = {\[Beta]r, \[Beta]i} /. 
    FindRoot[{Re[f1[\[Beta]r, \[Beta]i, \[Omega]a[[s]]]], 
      Im[f1[\[Beta]r, \[Beta]i, \[Omega]a[[s]]]]}, {{\[Beta]r, 
       Re[fa[[s]]]}, {\[Beta]i, Im[fa[[s]]]}}];
  \[Beta]a[[s]] = \[Beta]r1 + I*\[Beta]i1;
  ];
ListLinePlot[{\[Nu]a, Re@\[Beta]a}\[Transpose]]
ListLinePlot[{\[Nu]a, Im@\[Beta]a}\[Transpose]]

Вывод уравнения неоднократно проверялся.

[amon]

palz в сообщении #1421724 писал(а):

Вывод уравнения неоднократно проверялся.

А тут проверять особенно нечего. В цилиндрических координатах

$$\begin{align} E_z&=E_{z1}I_0(\varkappa_1r)\\ E_r&=\frac{ik_z}{\varkappa_1}E_{z1}I_1(\varkappa_1r)\quad\text{внутри провода}\\ E_z&=E_{z2}K_0(\varkappa_2r)\\ E_r&=\frac{ik_z}{\varkappa_2}E_{z2}K_1(\varkappa_2r)\quad\text{снаружи}\\ \varkappa_1&=\sqrt{k_z-\frac{\omega^2}{c^2}\varepsilon_1}\\ \varkappa_2&=\sqrt{k_z-\frac{\omega^2}{c^2}\varepsilon_2} \end{align}$$

Из граничных условий

$$\frac{\varepsilon_2}{\varkappa_2}I_0(\varkappa_1a)K_1(\varkappa_2a)-\frac{\varepsilon_1}{\varkappa_1}I_1(\varkappa_1a)K_0(\varkappa_2a)=0,$$

что знаком отличается от того, что Вы пытаетесь решить.

[palz]

Действительно, уравнение Гельмгольца приводит к $E_z(r)=E_{z1}I_0(\kappa_1 r)$ внутри и $E_z(r)=E_{z2}K_0(\kappa_2 r)$ снаружи провода. Далее ввиду отсутствия свободных зарядов можем написать $divD=0$. Поскольку $\varepsilon$ не зависит от координат, имеем и $divE=0$. В цилиндрических координатах это приводит к уравнению

$\frac{\partial E_r}{\partial r}+\frac{E_r}{r}+i\beta E_z=0$.

Здесь $E_\varphi=0$, т.к. рассматривается TM-мода. Решение данного уравнения ищем в виде $E_r=E_{r1}I_1(\kappa_1 r)$ внутри и $E_r=E_{r2}K_1(\kappa_2r)$ снаружи провода ($\kappa_1=\sqrt{\beta^2-\varepsilon_1k_0^2}$, $\kappa_2=\sqrt{\beta^2-k_0^2}$, $\beta$ - постоянная распространения, $k_0=\omega/c$). С учетом того, что

$ I_1'(z)=I_0(z)-\frac{1}{z}I_1(z), $
$ K_1'(z)=-K_0(z)-\frac{1}{z}K_1(z) $

Имеем

$ E_r=-\frac{i\beta}{\kappa_1}E_{z1}I_1(\kappa_1 r)$ (1) (здесь у вас нет минуса!)
$ E_z=E_{z1}I_0(\kappa_1 r)$, $r<a$ (2);
$ E_r=\frac{i\beta}{\kappa_2}E_{z2}K_1(\kappa_2 r) $ (3)
$E_z=E_{z2}K_0(\kappa_2 r)$, $r>a$ (4).

Граничные условия дают равенство компонент $E_z$ и $\varepsilon E_r$ внутри и снаружи провода. Учитывая, что снаружи $\varepsilon=1$, получим

$ \frac{\varepsilon_1 I_1(\kappa_1 a)}{\kappa_1 I_0(\kappa_1 a)}+\frac{K_1(\kappa_2 a)}{\kappa_2 K_0(\kappa_2 a)}=0 $