# Numerical scheme to calculate the normal mode of a set of hyperbolic PDEs?

I would like to solve the linearised, ideal, MHD equations, where the gas pressure is zero. $$\frac{\partial u_x}{\partial t}=v_A^2(x,z)\left[\nabla_{||}b_x - \frac{\partial b_{||}}{\partial x}\right],$$ $$\frac{\partial u_\perp}{\partial t}=v_A^2(x,z)\left[\nabla_{||}b_\perp - \nabla_\perp b_{||} \right],$$ $$\frac{\partial b_x}{\partial t}=\nabla_{||}u_x,$$ $$\frac{\partial b_\perp}{\partial t}=\nabla_{||}u_\perp,$$ $$\frac{\partial b_{||}}{\partial t}=-\left[\frac{\partial u_x}{\partial x}+\nabla_\perp u_\perp\right],$$ where $$\nabla_\perp = \cos\alpha\frac{\partial}{\partial y} - \sin\alpha\frac{\partial }{\partial z}$$ $$\nabla_{||} = \sin\alpha\frac{\partial}{\partial y} + \cos\alpha\frac{\partial }{\partial z}$$ and $$v_A(x,z)$$ denotes the background Alfven speed. The magnetic field perturbations $$b_x$$, $$b_\perp$$, $$b_{||}$$ are dimensionless because they have been normalised by $$B_0$$ and the background field is given by $$\boldsymbol{B}_0=B_0(\sin\alpha\,\boldsymbol{\hat{y}}+\cos\alpha\,\boldsymbol{\hat{z}}).$$

I would like to assume a $$y$$ and $$t$$ dependence of the form $$\exp[i(k_y y + \omega t)]$$ as this means we don't need to integrate in $$y$$ and $$t$$ which saves time.

Do you know a numerical scheme I could use to solve these equations? The boundary conditions I had in mind are periodic in $$z$$ and imposing a driver at $$x=x_{min}$$ and $$x=x_{max}$$, where $$v_A(x,z)$$, $$k_y$$ and $$\omega$$ are prescribed beforehand.

• What is the dimensionality of the problem? Are you trying to find just one eigenmode or several, or the full spectrum? There are multiple ways to approach this, the recommendations would depend on what numerical methods and computational software you are comfortable with. – Maxim Umansky Jun 18 '20 at 4:11
• I am not sure what you mean by dimensionality, the variables I would like to solve for are $u_x$, $u_\perp$, $b_x$, $b_\perp$ and $b_{||}$. These are all scalar functions of $x$ and $z$. I am actually wanting to prescribe $\omega$ and $k_y$ beforehand, impose a driver at $x=x_{min}$ and $x=x_{max}$ (periodic in $z$ say) then calculate the corresponding solution. A solution should exist provided I don't choose $k_y$ and $\omega$ in such a way that the solution contains singularities. – Peanutlex Jun 18 '20 at 9:44
• Ok, then write these equations discretized by finite-difference in Python so they are converted to a set of ODEs (method of lines), and use LSODA from scipy.integrate; that will work. I would recommend trying this first on a simpler but structurally similar problem, the acoustic wave in 1D, to build confidence – Maxim Umansky Jun 18 '20 at 14:58
• It will be stable with an implicit time-integrator like LSODA. You can use central differences. You can use low-order discretization, it does not matter too much at the initial stage. – Maxim Umansky Jun 19 '20 at 13:32
• LSODA as an adaptive implicit integrator would simplify requirements for your spatial discretization, so that would make it easier for initial development stage. Later on, once you understand better your problem and how the numerical solutions behave you could find that some other time-integration algorithm is more efficient here. But for getting something up and running quickly, without having much background in this area, LSODA is the way to go. It does not support complex numbers; but the way your equations are written there won't be any complex numbers there. – Maxim Umansky Jun 19 '20 at 20:50