首页
默认分类
技术经验
工作学习
娱乐爱好
闲言碎语
更多
统计
关于
登录
1
李芒果空岛-1.20.1-发展记录-05
191 阅读
2
“日晕“
183 阅读
3
初试3D打印——手机支架
139 阅读
4
使用typecho搭建博客
134 阅读
5
108第一届中国象棋比赛
129 阅读
Search
标签搜索
天文
Minecraft
李芒果空岛
数值计算
非线性最小二乘
typecho
macOS
PTCG
迭代法
Fortran
Halo
朗谬尔波
Langmiur
环法自行车赛
LM算法
博客框架
Hexo
Wordpress
3D打印
博客需求分析
Washy
累计撰写
62
篇文章
累计收到
1
条评论
首页
栏目
默认分类
技术经验
工作学习
娱乐爱好
闲言碎语
页面
统计
关于
管理后台
搜索到
8
篇与
的结果
2024-01-12
空间物理数据材料整理
0 前言 简单汇总下空间物理学领域部分数据、材料相关的FTP、网站等。 1 数据汇总 序号 简介 链接 备注 1 地磁与太阳活动指数 NOAA-FTP | GFZ-FTP | WDC 2 ACE卫星数据 NOAA-FTP 3 GIM Map GIPP-FTP 4 GOLD数据,地球同步轨道远紫外成像仪 GOLD 5 ESA Swarm卫星数据 Swarm 6 NASA太阳观测图像 NASA 7 SEPC太阳H$\alpha$图像 SEPC 8 测高仪数据 GIRO 2 材料汇总 序号 简介 网址 备注 1 国际电联通信相关文档 ITU-R 参考材料 地磁活动指数与太阳活动指数 GFZ数据下载的一种方式分享
2024年01月12日
49 阅读
0 评论
0 点赞
2023-09-25
史瓦西黑洞最内稳定圆轨道计算
0 前言 最内稳定圆轨道(Innermost Stable Circular Orbit,ISCO)是测试粒子可以稳定地绕广义相对论中的大质量物体运行的最小边缘稳定圆轨道。ISCO的半径$r_{\rm isco}$取决于中心物体的质量和角动量(自旋),它标记了黑洞吸积盘的内边缘。 ISCO 不应与罗希极限混淆,罗希极限是物理物体在潮汐力将其破坏之前可以绕轨道运行的最内点。ISCO 关注的是理论测试粒子,而不是真实物体。一般来说,ISCO 将比罗氏极限更接近中心物体。 本文介绍下如何推导史瓦西黑洞的最内稳定圆轨道,是一篇学习笔记。 1 有效势 在自然坐标系下($c=G=1$),史瓦西度规的球坐标形式如下 $$ d s^2 = -f dt^2 + \frac{d r^2}{f} + r^2 d \theta^2 + r^2 \sin^2 \theta d \phi \tag{1} $$ 其中$f = 1 - 2M/r$,$M$表示黑洞质量。因此度规张量形式为 $$ g_{\mu\nu} = \begin{bmatrix} -f & 0 & 0 & 0 \\ 0 & f^{-1} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2 \sin^2 \theta \end{bmatrix} $$ 考虑黑洞周围一个质量为$m$的测试粒子,其拉格朗日作用量为 $$ \widetilde{L} = \frac{m}{2} g_{\mu \nu} \dot{ x }^\mu \dot{ x}^\nu = \frac{m}{2} \left( -f \dot{t}^2 + f^{-1} \dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2 \theta \dot{\phi}^2 \right) $$ 为了计算方便,选取坐标系使得测试粒子的轨道面为赤道面,即$\theta = \pi/2$,则拉格朗日作用量简化为 $$ \widetilde{L} = \frac{m}{2} \left( -f \dot{t}^2 + f^{-1} \dot{r}^2 + r^2 \dot{\phi}^2 \right) \tag{2} $$ 根据四动量的定义$P_\mu = \partial \widetilde{L} / \partial \dot{ x}_\mu$,可知简化后的形式如下 $$ \begin{split} &P_t = \frac{\partial \widetilde{L}}{\partial \dot{t}} = -m f \dot{t} \equiv -E, \\ &P_r = \frac{\partial \widetilde{L}}{\partial \dot{r}} = mf^{-1} \dot{r}, \\ &P_\phi = \frac{\partial \widetilde{L}}{\partial \dot{\phi}} = m r^2 \dot{\phi} \equiv L. \end{split} \tag{3} $$ 其中$E$表示能量,$L$表示角动量,是两个守恒量。测试粒子的哈密顿量(Hamiltonian)为 $$ \begin{split} \widetilde{H} &= \dot{ x}^\mu P_\mu - \widetilde{L} = \frac{m}{2} g_{\mu \nu} \dot{ x}^\mu \dot{ x}^\nu \\ &= \frac{m}{2} \left( -f \dot{t}^2 + f^{-1} \dot{r}^2 + r^2 \dot{\phi}^2 \right) \end{split} \tag{4} $$ 由$\widetilde{H} = -mk/2$可知 $$ \dot{r}^2 = f(f\dot{t}^2 - r^2 \dot{\phi}^2 - k) = f \left( \frac{E^2}{m^2 f} - \frac{L^2}{m^2 r^2} - k \right) \tag{5} $$ 对上式进行整理,写为“动能”与“势能”之和的形式,即 $$ m^2 \dot{r}^2 + V(r) = E^2 \tag{6} $$ 其中$V(r)$为有效势,形式如下 $$ V(r) = f \left( \frac{L^2}{r^2} + m^2 k \right) \tag{7} $$ 2 圆轨道 测试粒子的运动是圆轨道的条件为 $$ \dot{r} = 0 \quad {\rm and} \quad \ddot{r} = 0 \tag{8} $$ 根据定义式$\ddot{r} \equiv d \dot{r} / dt$以及$\dot{r} \equiv dr/dt$可知 $$ \ddot{r} = \frac{d \dot{r}}{d t} = \frac{d \dot{r}}{dr} \frac{dr}{dt} = \dot{r} \frac{d \dot{r}}{dr} = \frac{1}{2} \frac{d \dot{r}^2}{dr} \tag{9} $$ 联立(5)式和(9)式可得 $$ \ddot{r} = 2ff'\dot{t}^2 - r(2f+rf') \dot{\phi}^2 - kf' \tag{10} $$ 将(5)式和(10)式代入(8)式可得 $$ \dot{t}^2 = \frac{2k}{2f - rf'}, \qquad \dot{\phi}^2 = \frac{kf'}{r(2f - rf')}. \tag{11} $$ 联立(3)式可知 $$ E^2 = \frac{2m^2kf^2}{2f - rf'}, \qquad L^2 = \frac{m^2 k r^3 f'}{2f-rf'}. \tag{12} $$ 3 最内稳定圆轨道 对于最内稳定圆轨道,除了圆轨道条件外,还需要满足有效势的二阶导数为零,即 $$ V''(r) = 0 \tag{13} $$ 首先,对$V(r)$求一阶导数,可得 $$ V'(r) = \frac{dV(r)}{dr} = f' \left( \frac{L^2}{r^2} + m^2k \right) + f \left( \frac{2LL'}{r^2} - \frac{2L^2}{r^3} \right) $$ 根据公式(3)的第三个公式可知$L' = 2L/r$,代入上式可得 $$ V'(r) = f' \left( \frac{L^2}{r^2} + m^2k \right) + \frac{2fL^2}{r^3} = \frac{2f + rf'}{r^3}L^2 + m^2kf' $$ 继续求导,有 $$ \begin{split} V''(r) &= -\frac{3(2f + rf')L^2}{r^4} + \frac{(2f' + f' + rf'')L^2 + (2f+rf')2LL'}{r^3} + m^2kf'' \\ &= (2f + 4rf' + r^2 f'') \frac{L^2}{r^4} + m^2kf'' \\ \end{split} \tag{14} $$ 将(12)式和(14)式代入(13)式,可得最内稳定圆轨道的大小为 $$ r_{\rm isco} = 6M $$ 参考 Innermost stable circular orbit https://ned.ipac.caltech.edu/level5/March01/Carroll3/Carroll7.html
2023年09月25日
56 阅读
0 评论
0 点赞
2023-09-14
高频传播特性的等效关系
0 前言 短波通信是一种常见的通信方式,其利用高频(High Frequency,HF)电磁波在电离层和地面之间的一次或多次反射,最远可将信号传至上万公里。在这个过程中,可以认为电离层和地面形成了类似“波导管”的结构。 在短波通信中,为了简化计算,通常会使用一些等效关系,这里将根据参考文献[1]中的内容,对这些等效关系进行介绍。 1 平面地球和平面电离层 1.1 “割线定律” - The "Secant Law" 在本节中,我们将考虑两个波的频率、虚拟路径和吸收之间存在的关系,一个波以斜入射反射,另一个波以垂直入射反射,二者的反射发生在相同的真实高度。为此,考虑以如下图所示的角度入射到平面电离层上的射线,其中电子密度随高度增加,从而发生全内反射。 Refs [1] 图4.1 平面地球与平面电离层的等效定理 在没有碰撞和外加磁场的情况下,在等离子体频率为$f_N$的水平高度上,频率为$f$的波的折射率$\mu$由下式给出: $$ \mu^2 = 1 - \left( \frac{f_N}{f} \right)^2 \tag{1} $$ 在反射层应用斯涅尔定律(Snell's law),即$\mu = \sin \phi_0$,可知 $$ f_N = f \cos \phi_0 $$ 如果$f_\nu$是与斜入射频率$f$在相同的真实高度(即相同的等离子体频率)上以垂直入射反射的频率,则$f_\nu =f_N$。因此 $$ f_\nu = f \cos \phi_0 \qquad {\rm or} \qquad f = f_\nu \sec \phi_0 \tag{2} $$ 该频率称为“等效垂直入射频率”,对应于$f$。上式中第二个式子就是所谓的割线定律。上述结果表明,在斜入射下,电离层反射的频率比正常入射时高得多。 1.2 Breit和Tuve定理 - Breit and Tuve's Theorem 另一个比较重要的关系称为Breit和Tuve定理。其中表示发射机T与接收机R之间传输的群(或等效)路径$P'$由等效三角形TAR的长度给出,即: $$ P' = TA + AR \tag{3} $$ 这可以通过以下论证来证明: $$ P' = \int_{TBR} \frac{ds}{\mu} = \int \frac{dx}{\mu \sin \phi} = \frac{1}{\sin \phi_0} \int dx = \frac{TR}{\sin \phi_0} = TA + AR $$ 请注意,反射的真实高度B总是小于A处的等效高度。 应该记住,只有当发送者和接收者位于电离层之外时(3)式才成立。如果发射机和接收机位于电离层内,即折射率为$\mu_1$的水平高度上,则(3)式的右侧必须除以$\mu_1$,这样$P'$仍然意味着群传播时间乘以自由空间速度。 1.3 马丁(等效路径)定理 - Martyn's (Equivalent Path) Theorem 如果$f$和$f_\nu$是从同一实际高度斜向和垂直反射的波的频率,则垂直反射信号的虚高等于斜向信号的等效三角路径的高度。 考虑等离子体频率为$f_N$的相同实际高度下,斜波和垂直波的折射率为$\mu_{ob}$和$\mu_\nu$,我们有 $$ \mu_{ob}^2 = 1 - \left( \frac{f_N}{f} \right)^2 \qquad {\rm and} \qquad \mu_\nu^2 = 1 - \left( \frac{f_N}{f \cos \phi_0} \right)^2 \tag{4} $$ 根据斯涅尔定律(Snell's law)可知$\mu_{ob} \sin \phi = \sin \phi_0$。联立这些方程可以得到 $$ \mu_{ob} \cos \phi = \mu_\nu \cos \phi_0 \tag{5} $$ 斜测信号的群路径为 $$ P' = \int_{TBR} \frac{ds}{\mu_{ob}} = 2 \int_0^{h_r} \frac{dh}{\mu_{ob} \cos \phi} = \frac{2}{\cos \phi_0} \int_0^{h_r} \frac{dh}{\mu_\nu} = \frac{2}{\cos \phi_0} h_\nu' $$ 其中 $$ h_\nu' = \frac{1}{2} P' \cos \phi_0 = \frac{1}{2} (TA + AR) \cos \phi_0 = AD $$ 因此 $$ P'(f) = h'(f_\nu) \sec \phi_0 \tag{6} $$ 这个定理表达了斜入射波的虚反射高度与等效垂直波的虚反射高度相等的重要关系。 1.4 马丁(吸收)定理 - Martyn's (Absorption) Theorem 略。 $$ \tag{7} $$ $$ \tag{8} $$ 2 电离层曲率的影响 Refs [1] 图4.2 斜入射射线几何图像 在弯曲的电离层中,斯涅耳定律的形式是 $$ \mu r \sin i = \mu_0 r_0 \sin i_0 \tag{9} $$ 式中$r$为从地心到折射率为$\mu$处的半径矢量的长度,半径矢量与光线夹角为$i$,如上图所示,图中$\mu_0$、$r_0$和$i_0$为任意参考值。以地面为参照,则有$\mu_0 = 1$、$r_0 = a$且$i_0 = (\pi/2) - \Delta$,其中$\Delta$为仰角。即 $$ \mu r \sin i = a \cos \Delta \tag{10} $$ 将式(1)代入式(10),使用$f_\nu$替换$f_N$,给出了频率$f$与相同实际高度$h_r$反射的等效垂直频率$f_\nu$之间的关系。即 $$ \left( \frac{f_\nu}{f} \right)^2 = 1 - \left( \frac{a \cos \Delta}{a + h_r} \right)^2 \tag{11} $$ 设$\phi_r'$为未折射光线的延线与半径矢量在$h_r$处的夹角,如上图所示。从几何关系中可以得到 $$ (a + h_r) \sin \phi_r = a \cos \Delta \tag{12} $$ 因此 $$ f_\nu = f \cos \phi_r \tag{13} $$ 由式(11)和式(13)可以看出,等效频率不仅与$\Delta$有关,而且与反射高度$h_r$有关。 $f_\nu$可以用(6)式或(7)式来定义,而不是用反射高度来定义,但与平面电离层的情况不同,$f_\nu$的值将取决于定义。 参考 [1] Ionospheric radio propagation
2023年09月14日
78 阅读
0 评论
0 点赞
2023-07-17
使用标势推导朗谬尔波
0 准备工作 引入标势$\varphi$和矢势$\mathbf{A}$,有 $$ \mathbf{E} = -\nabla \varphi, \qquad \mathbf{B} = \nabla \times \mathbf{A} $$ 则麦克斯韦方程组可以写为 $$ \nabla^2 \varphi = - \rho / \varepsilon_0, \qquad \nabla^2 \mathbf{A} - \nabla(\nabla \cdot \mathbf{A}) = - \mu_0 \mathbf{J} + \mu_0 \varepsilon_0 \nabla\frac{\partial \varphi}{\partial t} $$ 1 求解色散关系 1.1 双流程方程组 对于非磁化等离子体,双流体方程组的电子部分可以写为 $$ \left\{ \begin{split} &\nabla^2 \varphi = - \rho / \varepsilon_0 \\ &\frac{\partial n_e}{\partial t} + \nabla \cdot (n_e \mathbf{u_e}) = 0 \\ &n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = n_e e \nabla \varphi - \nabla p_e \end{split} \right. \tag{1} $$ 1.2 微扰法和平面波化 微扰法:离子作为正电荷背景,即$n_i = n_0$;设电子数密度由背景量和扰动量组成,即$n_e = n_0 + n_{e1}$。则电荷密度$\rho =-e n_{e1}$。 平面波化:设所有的扰动具有指数项$e^{i(\mathbf{k} \cdot \mathbf{r} - \omega t)}$,即具有平面波的形式,则有$\partial/\partial t = -i\omega$和$\nabla = i\mathbf{k}$。 因此双流体方程中电子部分可以改写为 $$ \left\{ \begin{split} &-k^2 \varphi = \frac{e}{\varepsilon_0} n_{e1} \\ &-i\omega n_{e1} + i n_0 \mathbf{k} \cdot \mathbf{u_e} = 0 \\ &-i\omega n_0 m_e \mathbf{u_e} = i \mathbf{k} n_0 e \varphi - i\mathbf{k} \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{2} $$ 运动速度$\mathbf{u_e}$只考虑波矢$\mathbf{k}$方向的分量。则上述方程组进一步简化为 $$ \left\{ \begin{split} &-k^2 \varphi = \frac{e}{\varepsilon_0} n_{e1} \\ &-\omega n_{e1} + k n_0 u_{e\parallel} = 0 \\ &-\omega n_0 m_e u_{e \parallel} = k n_0 e \varphi - k \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{3} $$ 联立上述方程组第二、第三两个方程,消去$u_{e\parallel}$可得 $$ n_{e1} = \frac{k^2 n_0 e}{k^2 \gamma_e k_B T_e - \omega^2 m_e} \varphi \tag{4} $$ 将上式代入方程组的第一个方程,消去$n_{e1}$和$\varphi$可得朗谬尔波的色散关系为 $$ \omega^2 = \frac{n_0 e^2}{\varepsilon_0 m_e} + \frac{\gamma_e k^2 k_B T_e}{m_e} = \omega_{pe}^2 + \frac{\gamma_e}{2} k^2 v_{the}^2 \tag{5} $$ 2 总结 从上述推导过程可以看出,使用标势推导色散关系相对而言更为简洁,且结论一致。标势的引入,也更有利于将结论推广至弯曲时空。此处,只是一个简单的尝试和复习以前的知识,并无独创之处。
2023年07月17日
57 阅读
0 评论
0 点赞
2023-06-25
朗谬尔振荡和朗谬尔波
1 朗谬尔(Langmiur)振荡 1.1 限制条件 考虑温度可以忽略不计的非磁化等离子体,对于与电子相关的现象(即高频部分),离子由于质量较大无法及时响应高频振荡,因此可看作不动的正电荷背景。当电子相对离子发生小扰动时,设扰动产生的扰动电场$\mathbf{E_1}$。 下图中红点表示离子,蓝点表示电子。 忽略磁场$\mathbf{B}$的等离子体称为非磁化等离子体,相对地,考虑磁场$\mathbf{B}$的等离子体称为磁化等离子体。 设$n_{e,i}$表示电子/离子数密度,$\mathbf{u_{e,i}}$表示电子/离子速度,$m_{e,i}$表示电子/离子质量,则电荷密度$\rho = e (n_i - n_e)$。 1.2 双流体方程组 由高斯定律可知 $$ \nabla \cdot \mathbf{E_1} = \frac{\rho}{\varepsilon_0} = \frac{e}{\varepsilon_0}(n_i - n_e) \tag{1} $$ 由法拉第电磁感应定律可知 $$ \nabla \times \mathbf{E_1} = - \frac{\partial \mathbf{B}}{\partial t} = 0 \tag{2} $$ 电子和离子的运动需要分开讨论,由于离子作为不动的正电荷背景,此处只考虑电子。由连续性方程可知 $$ \frac{\partial n_e}{\partial t} + \nabla \cdot (n_e \mathbf{u_e}) = 0 \tag{3} $$ 由运动方程可知 $$ n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = -n_e e(\mathbf{E_1} + \mathbf{u_e} \times \mathbf{B}) = - n_e e \mathbf{E_1} \tag{4} $$ 1.3 微扰法和平面波化 微扰法:离子作为正电荷背景,即$n_i = n_0$;设电子数密度由背景量和扰动量组成,即$n_e = n_0 + n_{e1}$。代入上述方程组,并忽略二阶小项,则方程组可改写为 $$ \left\{ \begin{split} &\nabla \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &\nabla \times \mathbf{E_1} = 0 \\ &\frac{\partial n_{e1}}{\partial t} + n_0 \nabla \cdot \mathbf{u_e} = 0 \\ &m_e \frac{\partial \mathbf{u_e}}{\partial t} = -e\mathbf{E_1} \end{split} \right. \tag{5} $$ 平面波化:设所有的扰动具有指数项$e^{i(\mathbf{k} \cdot \mathbf{r} - \omega t)}$,即具有平面波的形式。利用$\partial/\partial t = -i\omega$和$\nabla = i\mathbf{k}$对上述方程组进行化简 $$ \left\{ \begin{split} &i\mathbf{k} \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\mathbf{k} \times \mathbf{E_1} = 0 \\ &-i\omega n_{e1} + i n_0 \mathbf{k} \cdot \mathbf{u_e} = 0 \\ &-i\omega m_e \mathbf{u_e} = -e\mathbf{E_1} \end{split} \right. \tag{6} $$ 由$i\mathbf{k} \times \mathbf{E_1} = 0$可知$\mathbf{k} \parallel \mathbf{E_1}$,运动速度$\mathbf{u_e}$只考虑扰动电场$\mathbf{E_1}$方向的分量。则上述方程组进一步简化为 $$ \left\{ \begin{split} &i k E_1 = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\omega n_{e1} - i n_0 k u_{e\parallel} = 0 \\ &i\omega m_e u_{e\parallel} = e E_{1} \end{split} \right. \tag{7} $$ 由上述方程组的第二和第三式消去$u_{e\parallel}$可得 $$ n_{e1} = - \frac{i k n_0 e}{\omega^2 m_e} \tag{8} $$ 将(8)式代入(7)式第一个方程,消去$E_1$、$n_{e1}$可得 $$ \omega^2 = \frac{n_0 e^2}{\varepsilon_0 m_e} \Rightarrow \omega_{pe} = \sqrt{\frac{n_0 e^2}{\varepsilon_0 m_e}} \tag{9} $$ 其中$\omega_{pe}$表示朗谬尔频率,也叫做(电子)等离子体频率。由群速度$v_g = d \omega/dk = 0$可知,该振荡只存在于振荡产生的位置,不会向外传播。 2 朗谬尔波 2.1 限制条件 当电子温度不为零且其他条件不变时,由于电子的热运动,可以将振荡区域的信息携带至邻近区域,从而使邻近区域也发生振荡。这样,发生在某处的振荡就能传播出去而形成波,这种波称为等离子体波或朗谬尔波$^{[1]}$。 [1] 《等离子体物理学》李定著,P92. 2.2 色散关系 此时,高斯定律、法拉点电磁感应定律和连续性方程的形式不发生改变,运动方程中需要引入热压梯度项$\nabla p_e$,如下 $$ n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = -n_e e(\mathbf{E_1} + \mathbf{u_e} \times \mathbf{B}) - \nabla p_e = -n_e e \mathbf{E_1} - \nabla p_e \tag{10} $$ 由状态方程$p_e \rho_e^{-\gamma_e} = {\rm consts}$和$p_e = n_e k_B T_e$以及$\rho_e = -e n_e$,可得 $$ \begin{split} &\nabla(p_e \rho_e^{-\gamma_e}) = \rho_e^{-\gamma_e} \nabla p_e + p_e \nabla \rho_e^{-\gamma_e} = \rho_e^{-\gamma_e} \nabla p_e - p_e \gamma_e \rho_e^{-\gamma_e - 1} \nabla \rho_e = 0 \\ \Rightarrow & \nabla p_e = p_e \gamma_e \rho_e^{-1} \nabla \rho_e = \gamma_e k_B T_e \nabla n_{e1} \end{split} \tag{11} $$ 其中$k_B$表示玻尔兹曼常数,$T_e$表示电子温度,$\gamma_e$表示电子比热比(在绝热假设中,$\gamma_e = 3$;在等温假设中,$\gamma_e = 1$)。将上式代入(10)式,则完整的方程组可以写为 $$ \left\{ \begin{split} &\nabla \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &\nabla \times \mathbf{E_1} = 0 \\ &\frac{\partial n_{e1}}{\partial t} + n_0 \nabla \cdot \mathbf{u_e} = 0 \\ &n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = -n_e e \mathbf{E_1} - \gamma_e k_B T_e \nabla n_{e1} \end{split} \right. \tag{12} $$ 平面波化后,上述方程组改写为 $$ \left\{ \begin{split} &i\mathbf{k} \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\mathbf{k} \times \mathbf{E_1} = 0 \\ &-i\omega n_{e1} + i n_0 \mathbf{k} \cdot \mathbf{u_e} = 0 \\ &-i\omega n_0 m_e \mathbf{u_e} = -n_0 e\mathbf{E_1} - i k \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{13} $$ 同样地,运动速度$\mathbf{u_e}$只考虑扰动电场$\mathbf{E_1}$方向的分量,将上述方程组化简为 $$ \left\{ \begin{split} &i k E_1 = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\omega n_{e1} - i n_0 k u_{e\parallel} = 0 \\ &i\omega n_0 m_e u_{e\parallel} = n_0 e E_{1} + i k \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{14} $$ 由上述方程组的第二和第三式消去$u_{e\parallel}$可得 $$ n_{e1} = - \frac{ik n_0 e}{\omega^2 m_e - k^2 \gamma_e k_B T_e} E_1 \tag{15} $$ 将(15)式代入(14)式第一个方程,消去$E_1$、$n_{e1}$可得 $$ \omega^2 = \frac{n_0 e^2}{\varepsilon_0 m_e} + \frac{\gamma_e k^2 k_B T_e}{m_e} = \omega_{pe}^2 + \frac{\gamma_e}{2} k^2 v_{the}^2 \tag{16} $$ 其中$v_{the} = \sqrt{2 k_B T_e / m_e}$表示电子热速度。上式即为朗谬尔波的色散关系,从(16)式可以看出,只有当$\omega > \omega_{pe}$时,朗谬尔波才能传播。
2023年06月25日
89 阅读
0 评论
0 点赞
2023-03-29
非线性最小二乘法拟合函数-3
在前面两个博客中,推导得到了Levenberg-Marquardt算法(简称LM算法)的迭代公式,这里将讲述如何使用Fortran编写一个简单的LM算法。 此处介绍待定系数只有一个的情况。 1 程序设计 1.1 Levenberg_Maquardt_Fit 简介:LM算法子例程,输入离散数据,迭代拟合待定系数 传递参数: n 离散数据的长度 x 离散数据的自变量 y 离散数据的因变量 a 待定系数 1.2 myfunc 简介:待拟合函数的数值计算子例程。为了具有普适性,假定待拟合函数的解析式未知,函数结果只能由此子例程数值计算得到。 传递参数: n 离散数据的长度 x 离散数据的自变量 a 待定系数的当前值 fx 待拟合函数的数值计算结果 1.3 Calculate_Jacobian 简介:计算待拟合函数的Jacobian矩阵的子例程。为了具有普适性,使用差分法计算一阶导数。 传递参数: n 离散数据的长度 x 离散数据的自变量 y 待拟合函数的数值计算结果 a 待定系数的当前值 J 待拟合函数的Jacobian矩阵 2 源代码 2.1 Levenberg_Maquardt_Fit 此处迭代次数上限设置为30次,可根据需求自行更改 对于自变量有多个的情况,可将输入x改为x1,x2,...,满足自变量个数即可 SUBROUTINE Levenberg_Maquardt_Fit(n, x, y, a) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n ! 数据点个数 REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点 ! REAL, DIMENSION(3), INTENT(INOUT) :: a ! 待求的拟合系数 REAL, INTENT(INOUT) :: a ! 待求的拟合系数 ! 定义常量 INTEGER, PARAMETER :: m = 1 ! 待求的系数个数 REAL, PARAMETER :: eps = 1.0E-6 ! 收敛阈值 ! 定义变量 REAL :: da(m), r(n), J(n,m), H(m,m) ! 拟合系数、残差、雅可比矩阵、Hessian矩阵 REAL :: lambda, alpha(m) ! 调节因子、步长 INTEGER :: i, iter ! 循环计数器 REAL :: fx(n) ! 初始化调节因子 lambda = 0.001 ! 开始迭代 iter = 0 DO WHILE(iter < 30) ! 迭代次数上限为10000 iter = iter + 1 ! 计算被拟合函数值 CALL myfunc(n,x,a,fx) ! 计算残差向量 r = y - fx ! 计算雅可比矩阵 CALL Calculate_Jacobian(n, x, fx, a, J) ! 计算Hessian矩阵 H = MATMUL(TRANSPOSE(J), J) ! 计算梯度向量 da = MATMUL(TRANSPOSE(J), r) ! 计算搜索方向 DO i = 1, m H(i,i) = H(i,i) + lambda END DO ! CALL SOLVE_LINEAR_SYSTEM(m, H, da, alpha) alpha = da(1)/H(1,1) ! 计算新的拟合系数 a = a + alpha(1) ! 更新调节因子 IF (NORM2(r) < eps) THEN EXIT ELSE IF (NORM2(r) < NORM2(r - MATMUL(J, alpha))/2) THEN lambda = lambda/10. ELSE lambda = lambda*10. END IF END DO END SUBROUTINE 2.2 myfunc 此处使用了一个非常简单的一次函数作为示例,使用时需要将fx = a*x修改为实际的待拟合函数的形式。若待拟合函数是一个子例程,则在此处使用CALL调用 若在Levenberg_Maquardt_Fit子例程中修改了自变量个数,此处需要对应修改 SUBROUTINE myfunc(n,x,a,fx) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x(n) ! 自变量 REAL, INTENT(IN) :: a ! 拟合系数 REAL, INTENT(OUT) :: fx(n) fx = a*x END SUBROUTINE 2.3 Calculate_Jacobian 本子例程使用了差分法计算一阶导数,因此几乎适用于所有的待拟合函数 若在前面修改了自变量的个数,此处需对应修改 SUBROUTINE Calculate_Jacobian(n, x, y, a, J) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n ! 数据点个数 REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点 REAL, INTENT(IN) :: a ! 拟合系数 ! 输出参数 REAL, DIMENSION(n,1), INTENT(OUT) :: J ! 雅可比矩阵 ! 定义变量 INTEGER :: i ! 循环计数器 REAL :: dela, y1(n) ! 步进 dela = 0.01 ! 计算 a+dela 对应的 y1 CALL myfunc(n,x,a+dela,y1) ! 计算雅可比 DO i = 1, n J(i,1) = (y1(i) - y(i))/dela END DO END SUBROUTINE 3 测试 编写程序主体,设定y约为x的两倍,调用LM子例程计算待定系数a program name implicit none INTEGER n real :: x(5),y(5),a x = (/1.,2.,3.,4.,5./) y = (/2.1,4.1,5.9,8.1,9.9/) ! y = 2*x**2 + n = 5 a = 1.0 call Levenberg_Maquardt_Fit(n, x, y, a) print *,a end program name 使用gfortran编译生成可执行文件,运行可执行文件后,在终端打印出 1.99818182 离散数据的因变量约为自变量的两倍,因此上述计算结果符合预期
2023年03月29日
72 阅读
0 评论
0 点赞
2023-03-29
非线性最小二乘法拟合函数-2
问题1:假设有$n$个离散点$(x_i,y_i)$,存在一个系数未知的被拟合函数$f(x,a)$,如何使用非线性最小二乘法拟合得到最优的待定系数$a$。 问题2:假设有$n$个离散点$(x_i,y_i)$,存在一个系数未知的被拟合函数$f(x,a_1,a_2,\dots,a_m)$,如何使用非线性最小二乘法拟合得到最优的待定系数组$a_1,a_2,\dots,a_m$。 1 数学推导 上一个博客讲解了问题1的数学推导,这里关注待定系数有多个的问题2。 1.1 问题数学化 对于问题2,第$i$个数据的残差表达式如下 $$ r_i (a_1,a_2,\dots,a_m) = y_i - f(x_i,a_1,a_2,\dots,a_m) = y_i - f_i(a_1,a_2,\dots,a_m) $$ 上式中,$x_i,y_i$均为已知量,为书写方便,将$f(x_i,a_1,a_2,\dots,a_m)$简写为$f_i(a_1,a_2,\dots,a_m)$。 使用残差的平方和作为损失函数,如下 $$ L(a_1,a_2,\dots,a_m) = \sum_{i=1}^n r_i^2(a_1,a_2,\dots,a_m) = \sum_{i=1}^n[y_i - f_i(a_1,a_2,\dots,a_m)]^2 $$ 则问题2可表示为求解损失函数$L(a_1,a_2,\dots,a_m)$最小值对应的待定系数组$a_1,a_2,\dots,a_m$。同样的,将最小值问题近似为极小值问题。为书写方便,使用向量$\mathbf{a}$表示$a_1,a_2,\dots,a_m$,则上式修改为 $$ L(\mathbf{a}) = \sum_{i=1}^n [y_i - f_i(\mathbf{a})]^2 $$ 1.2 求函数极值点 1.2.1 迭代法递推公式 求解函数的极值点问题,即为求解一阶导数的零点问题,求函数$L(\mathbf{a})$对$\mathbf{a}$的一阶导数,如下 $$ \frac{d L(\mathbf{a})}{d \mathbf{a}} = -2 \sum_{i=1}^n [y_i - f_i(\mathbf{a})] \frac{d f_i(\mathbf{a})}{d \mathbf{a}} = -2 \sum_{i=1}^n [y_i - f_i(\mathbf{a})] J_i(\mathbf{a}) $$ 其中$J_i(\mathbf{a}) = d f_i(\mathbf{a}) / d \mathbf{a}$。令 $$ g(\mathbf{a}) = -\frac{1}{2} \frac{d L(\mathbf{a})}{d \mathbf{a}} = \sum_{i=1}^n J_i(\mathbf{a}) [y_i - f_i(\mathbf{a})] $$ 至此,将求解函数$L(\mathbf{a})$的极值点问题转化为求解函数$g(\mathbf{a})$的零点问题。对$g(\mathbf{a})$在$\mathbf{a}=\mathbf{a_0}$处进行泰勒展开,其中$\mathbf{a_0}$表示待定系数组$\mathbf{a}$的初始值,只保留一阶项有 $$ g(\mathbf{a}) \approx \sum_{i=1}^n J_i(\mathbf{a_0}) [y_i - f_i(\mathbf{a_0}) - J_i(\mathbf{a_0}) (\mathbf{a} - \mathbf{a_0})] = \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) - \mathbf{H} (\mathbf{a_0}) \Delta \mathbf{a} $$ 其中$\mathbf{H} (\mathbf{a}) = \mathbf{J^T} \mathbf{J}$为Hession矩阵,$\mathbf{J}(\mathbf{a})$为Jacobian矩阵,表达式如下 $$ \mathbf{J} (\mathbf{a}) = \left[ \begin{matrix} \frac{\partial J_1}{\partial a_1} & \frac{\partial J_1}{\partial a_2} & \cdots & \frac{\partial J_1}{\partial a_m} \\ \frac{\partial J_2}{\partial a_1} & \frac{\partial J_2}{\partial a_2} & \cdots & \frac{\partial J_2}{\partial a_m} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial J_n}{\partial a_1} & \frac{\partial J_n}{\partial a_2} & \cdots & \frac{\partial J_n}{\partial a_m} \end{matrix} \right] $$ 令$g(\mathbf{a})=0$可得 $$ \mathbf{H} (\mathbf{a_0}) \Delta \mathbf{a} = \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 上式即为高斯-牛顿法迭代公式。引入阻尼因子$\mu$得到Levenberg-Marquardt算法迭代公式如下 $$ [ \mathbf{H} (\mathbf{a_0}) + \mu \mathbf{I} ] \Delta \mathbf{a} = \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 其中$\mathbf{I}$为单位矩阵。整理后可得 $$ \Delta \mathbf{a} = [ \mathbf{H} (\mathbf{a_0}) + \mu \mathbf{I} ]^{-1} \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 进而可以得到递推公式为 $$ \mathbf{a_{k+1}} = \mathbf{a_{k}} + [ \mathbf{H} (\mathbf{a_0}) + \mu \mathbf{I} ]^{-1} \mathbf{J^T}(\mathbf{a_0}) \mathbf{r} (\mathbf{a_0}) $$ 至此,得到了非线性最小二乘拟合的Levenberg-Marquardt算法迭代公式。 上述迭代公式得到的是待定系数组。 1.2.2 Jocabian矩阵 对于可直接求导得到一阶导数解析表达式的函数,直接代入相应的解析表达式即可。对于只能数值求解的函数,可使用差分近似求解,即 $$ J(\mathbf{a}) \approx \frac{f(\mathbf{a}+\delta \mathbf{a}) - f(\mathbf{a})}{\delta \mathbf{a}} $$ 其中$\delta \mathbf{a}$表示一个小量。至此,即可使用Levenberg-Marquardt算法迭代公式进行迭代求解。 1.2.3 结束迭代的条件 残差的范数小于某个值 超出最大循环次数 2 两个问题的联系 在本文中,Jacobian矩阵的大小为n行m列,Hession矩阵的大小为m行m列。其中,n为离散点的个数,m为待定系数的个数。当$m=1$时,问题2退化为问题1,即Jacobian矩阵变为n行1列的向量,Hession矩阵变为一个标量。
2023年03月29日
72 阅读
0 评论
0 点赞
2023-03-29
非线性最小二乘法拟合函数-1
从数学推导至程序实现,对非线性最小二乘拟合原理及应用进行简单的介绍,希望能够帮助到有需要的人。 0 问题引入 问题1:假设有$n$个离散点$(x_i,y_i)$,存在一个系数未知的被拟合函数$f(x,a)$,如何使用非线性最小二乘法拟合得到最优的待定系数$a$。 问题2:假设有$n$个离散点$(x_i,y_i)$,存在一个系数未知的被拟合函数$f(x,a_1,a_2,\dots,a_m)$,如何使用非线性最小二乘法拟合得到最优的待定系数组$a_1,a_2,\dots,a_m$。 问题2是问题1的拓展,由一个待定系数变为多个待定系数。之后的介绍和推导将由浅入深,从问题1开始讲起。 1 非线性最小二乘法 概念:非线性最小二乘法是一种优化方法,它通过最小化非线性函数与实验数据之间的残差平方和来确定未知参数的最优解。 非线性函数:参数未知的非线性被拟合函数,常用的有多项式函数、指数函数、幂函数、高斯函数等。 残差:是指在拟合过程中,模型预测值与实际观测值之间的差异。 对残差的要求: 独立同分布:残差之间应该是独立的且来自同一个概率分布。一般指高斯分布,如果不是,可能会出现偏差、误差较大的问题; 零均值:残差的平均值应该为0; 常数方差:残差的方差应该为常数。 求解方法:通常使用迭代优化算法来寻找最优解,如高斯-牛顿法、Levenberg-Marquardt算法等。 2 迭代法 概念:迭代法也称辗转法,是一种不断用变量的旧值递推新值的过程。迭代法是用计算机解决问题的一种基本方法,它利用计算机运算速度快、适合做重复性操作的特点,让计算机对一组指令(或一定步骤)进行重复执行,在每次执行这组指令(或这些步骤)时,都从变量的原值推出它的一个新值,迭代法又分为精确迭代和近似迭代。 百度百科——迭代法 下面介绍两种使用迭代法的场景。 2.1 求函数的零点 对函数$f(x)$在$x=x_0$处进行泰勒展开,只保留到一阶项,即 $$ f(x) \approx f(x_0) + f'(x_0) (x - x_0) $$ 令上式等于零,若函数$f(x)$的一阶导数存在且不为零,可得 $$ x = x_0 - f(x_0)/f'(x_0) $$ 进而可以得到递推公式 $$ x_{k+1} = x_k - f(x_k)/f'(x_k) $$ 设增量$\Delta = x_{k+1} - x_k$,则 $$ \Delta = - f(x_k)/f'(x_k) $$ 根据递推公式即可完成迭代求解,本方法为牛顿迭代法。 2.2 求函数的极值点 已知极值点处的一阶导数为零,即问题可修改为求函数$f'(x)$的零点。利用2.1节的方法,对$f'(x)$进行泰勒展开,可得 $$ f'(x) \approx f'(x_0) + f''(x_0) (x - x_0) $$ 同样的,可得增量表达式为 $$ \Delta = - f'(x_k)/f''(x_k) $$ 3 数学推导 由浅入深,先讲解待定系数只有一个的问题1。 3.1 问题数学化 问题1:假设有$n$个离散点$(x_i,y_i)$,存在一个系数未知的被拟合函数$f(x,a)$,如何使用非线性最小二乘法拟合得到最优的待定系数$a$。 对于问题1,第$i$个数据的残差表达式如下 $$ r_i (a) = y_i - f(x_i,a) = y_i - f_i(a) $$ 上式中,$x_i,y_i$均为已知量,为书写方便,将$f(x_i,a)$简写为$f_i(a)$。 根据第1节的概念,使用残差的平方和作为损失函数,表达式如下 $$ L(a) = \sum_{i=1}^n r_i^2(a) = \sum_{i=1}^n[y_i - f_i(a)]^2 $$ 则问题1可表示为求解损失函数$L(a)$最小值对应的待定系数$a$。最小值的问题通常不容易求解,因此可将其近似为求解损失函数$L(a)$极小值对应的待定系数$a$。若最终求得的极小值不是最小值,则为局部最优解。 3.2 求函数L(a)的极值点 3.2.1 迭代法递推公式 求解函数的极值点问题,即为求解一阶导数的零点问题,求函数$L(a)$对$a$的一阶导数,如下 $$ \frac{d L(a)}{d a} = -2 \sum_{i=1}^n [y_i - f_i(a)] \frac{d f_i(a)}{d a} $$ 令$J_i (a) = d f_i (a) / d a$,则有 $$ g(a) = -\frac{1}{2} \frac{d L(a)}{d a} = \sum_{i=1}^n [y_i - f_i(a)] J_i (a) $$ 至此,将求解函数$L(a)$的极值点问题转化为求解函数$g(a)$的零点问题。对$g(a)$在$a=a_0$处进行泰勒展开,并且只保留一阶项,有 $$ g(a) \approx \sum_{i=1}^n J_i (a_0) [y_i - f_i(a_0) - J_i(a_0) (a - a_0)] = \sum_{i=1}^n J_i (a_0) r_i(a_0) - H(a_0) \Delta a $$ 其中$H(a_0) = \sum_{i=1}^n J_i^2(a_0)$,$\Delta a = a - a_0$。令$g(a)=0$可得 $$ H(a_0) \Delta a = \sum_{i=1}^n J_i (a_0) r_i(a_0) $$ 上式即为高斯-牛顿法迭代公式。引入阻尼因子$\mu$得到Levenberg-Marquardt算法迭代公式如下 $$ [H(a_0) + \mu] \Delta a = \sum_{i=1}^n J_i (a_0) r_i(a_0) $$ 整理后可得 $$ \Delta a = [H(a_0) + \mu]^{-1} \sum_{i=1}^n J_i(a_0) r_i(a_0) $$ 进而可以得到递推公式为 $$ a_{k+1} = a_k + [H(a_k) + \mu]^{-1} \sum_{i=1}^n J_i(a_k) r_i(a_k) $$ 至此,得到了非线性最小二乘拟合的Levenberg-Marquardt算法迭代公式。 阻尼因子:是用来平衡算法的收敛速度和稳定性的一个参数。 当阻尼因子较小时,迭代算法会更快地收敛,但也会更容易陷入局部最小值; 当阻尼因子较大时,算法会更稳定,但会收敛得更慢。 3.2.2 一阶导数 对于可直接求导得到一阶导数解析表达式的函数,直接代入相应的解析表达式即可。对于只能数值求解的函数,可使用差分近似求解,即 $$ J(a) \approx \frac{f(a+\delta a) - f(a)}{\delta a} $$ 其中$\delta a$表示一个小量。至此,即可使用Levenberg-Marquardt算法迭代公式进行迭代求解。 3.2.3 结束迭代的条件 残差的范数小于某个值 超出最大循环次数
2023年03月29日
93 阅读
0 评论
0 点赞