Let $V$ be the vector space of trigonometric polynomials of degree $\le D$ on the flat torus $\mathbb T^n$. That is, $$V=\operatorname{Span}\left\{\cos (2\pi \lambda \cdot x), \sin (2\pi \lambda \cdot x)\mid\lambda\in\mathbb Z^n,|\lambda|_1 \le D\right\}$$
Are the Morse functions (functions for which zero is a regular value of the gradient) dense in $V$?
What would be the simplest argument for this (I know calculus, but I don't know Morse theory)?
Lets observe the $1$-dimensional case for a start. You can approximate cosines, just take $\cos(x)+\frac{1}{N}\sin(x)$ as a function with non-zero derivative approximating cosine by error bounded $O(1/N)$, and clearly all the sine functions have non-zero derivatives at $0$.
Now notice that if $f,g$ are two Morse-like functions, then their sum can be approximated as well in the following manner - if their sum has non-zero derivative, then the claim is obvious, but if the sum is $0$, then look at $(1+\frac{1}{n})f(x)+(1-\frac{1}{n})g(x)$, the difference from the original sum is bounded by $\frac{2\max\{\|f\|_{\infty},\|g\|_{\infty}\}}{N}$ in the sup. norm, and its derivative equals to $\frac{f'(0)-g'(0)}{N}\neq 0$.
Therefore, given a general polynomial expression, just approximate the basis functions and then approximate the whole sum (as your space is finite dimensional).
P.S. by changing the approximation rate with each summand, and taking rapidly decreasing sequence, one can approximate like that a Fourier series, say assuming you have nice convergence by assumption on the decay rate of the coefficients (and even keep the approximation in a suitable Sobolev space).