Total Harmonic Distortion calculation and its origins
You can create non-linear digital systems (an example would be a system that finds the absolute value of the input). You can also simulate an analog non-linear system using DSP. The easiest way is to use a power series. Assuming that the non-linear system you're modeling is time-invariant and memoryless, you can approximate it as follows:

y(t) = a0 + a1*x(t) + a2*x(t).*x(t) + a3*x(t).*x(t).*x(t) + ....


where y(t) is the system output, x(t) is the input, and an are coefficients that depend on the system you're modeling.

As an example, you can do a simple simulation of mild harmonic distortion by setting x(t) equal to the a cosine signal, and setting a1=2, a3=-0.1, and all other coefficients equal to 0. a1 is the system gain, and a3 is the non-linear gain. By setting a3 to larger negative values, you increase the nonlinear effect -- you can think of it as an amplifier getting closer to saturation.

If you plot the spectra of x(t) and y(t), you'll clearly see harmonic distortion occuring.