# Suitable finite difference method for a convection-diffusion system?

I am trying to solve a system of PDEs

$H_{t} = \frac{0.3}{0.7} - \frac{0.005 B f(h(H))}{\theta} - \frac{0.3 f(h(H))}{0.7} + \frac{500}{0.7} (HH_x)_x + (HH_y)_y$

$N_t = \frac{N_{in} - 0.002 [N] B f(h(H)) + 0.1 \frac{0.002}{g} B - 0.1 N}{0.7 H} + 10 ([N]_{xx} + [N]_{yy}) + \frac{500}{0.7} ([N] H_x)_x + ([N] H_y)_y$

$B_t = g [N] B f(h(H)) - 0.3 B + 2 (B_{xx} + B_{yy})$

using the method of lines. What finite difference method(s) would be suitable for attempting a solution? Most of the examples I can find are either in 1 spatial dimension only or for a single variable. Since this is a convection-diffusion system

• It would be easier if the equations looked more alike to convection-diffusion equations. – nicoguaro Nov 14 '17 at 16:20
• The equations have a nice divergence structure which may indicate some physical conservation principle. It may be better to keep this form and directly make a scheme for it. – cfdlab Nov 15 '17 at 3:53

The $B$ equation is parabolic and you can use standard central finite difference. The $H$ equation also looks like parabolic so you can do a central difference like $$(H H_x)_x = \frac{H_{i+1/2}(H_{i+1}-H_i) - H_{i-1/2}(H_i - H_{i-1})}{h^2}, \qquad H_{i+1/2}=\frac{H_i + H_{i+1}}{2}$$ Equivalently we could approximate as $$(HH_x)_x = \frac{(H^2)_{xx}}{2} = \frac{H_{i-1}^2 - 2H_i^2 + H_{i+1}^2}{2h^2}$$ If $H$ becomes zero or negative, that could be problematic since it corresponds to vanishing viscosity or negative viscosity. Does this happen in your model?
Is there no time derivative in the middle equation? For $N_{xx}$ etc., you can do central difference and terms like $(NH_x)_x$ can be again central differenced as in the first formula above. Note that these suggestions are just based on consistency and you have to look at the stability of these schemes.