Wikipedia's page on Axial Precession has a good deal of mathematics on it and, unless you have absolutely got to have precise numbers for some obscure story reason (which is probably what I'd call the tail wagging the dog - change the story to avoid that problem) then making up the numbers or just forgetting all about Axial Precession is the way to go. The periods involved will be very long by the standards of any story you're likely to write, so why burden yourself with something you probably don;t need.
The maths that follows only gives a rough approximation anyway, and you don't want to even contemplate the kind of things you need to do to get a better one : it's not worth it.
That said let's have a look at the very basic theory result that Wikipedia gives :
There are two component to axial precession that matter (for Earth) : the one due to the Moon and the one due to the Sun. The Moon's effect is actually larger, but these numbers are very sensitive to the values you use.
The Solar Contribution
$$\frac{d\psi}{dt}=\left[ \frac {GM_s}{a_s^3\left(1-e_s^2\right)^\frac 3 2} \right] \left[ \frac {C-A}C \frac {cos\epsilon}\omega \right]$$
Lots of symbols so what do they mean ?
- $G$ - the Universal Gravitational Constant also famous from $F = \frac {GM_1M_2} {r^2}$ Newtons law for gravitation.
- $M_s$ - The Sun's mass - in your case you need the mass of your planet's star, of course.
- $a_s$ - The semi-major axis of the orbit of the planet around it's star
- $e_s$ - The eccentricity of the planet's orbit around it's star.
Now that second term in square bracket, which is also in the expression for the Lunar contribution. This one is trickier.
- $C$ - moment of inertia (of Earth) around the axis of rotation
- $A$ - moment of inertia around the equator
- $\epsilon$ - the angle between the equatorial plane and the ecliptic plane (see below)
- $\omega$ - Earth's angular velocity (due to it's rotation, not it's orbit)
Now this expression is really poorly dealt with in Wikipedia because of two problems.
$epsilon$ in Wikipedia is assumed to be the same for both the Solar and Lunar contributions. This is not (AFAIK) correct. The angle should be the angle between axis of rotation of the body and the plane of the orbit of the other body (which means it's different for the Sun and Moon).
The $C-A$ and $A$ terms are really hard to deal with for mere mortals (and frankly just messy for anyone else). For your purposes I would propose the following compromise term instead. It's based on modeling the planetary bulge as an ellipsoid of constant density compared with the $A$ value for an ideal sphere - both objects have the same mass and density, which I'm taking as constant . I'll spare you the derivation :
$$\frac {C-A} C \approx 1 - \frac{R^2} {a^2}$$
where in this case $R$ is the average radius of the planet and $a$ is the equatorial radius of the planet.
The Lunar Contribution
$$\frac{d\psi}{dt}=\left[ \frac {GM_l}{a_l^3\left(1-e_l^2\right)^\frac 3 2} \left(1-\frac 3 2 sin^2i\right) \right] \left[ \frac {C-A}C \frac {cos\epsilon}\omega \right]$$
Not much change here expect that the masses and so on refer to the Moon and not the Sun (hence the different subscripts). There is one additional term which is the factor :
$$\left(1-\frac 3 2 sin^2i\right)$$
This corrects for the effect that the angle of inclination of Moon's orbit to the ecliptic is not zero. The ecliptic being the plane with the Sun and the Earth's orbit in it.
You have to decide these numbers for yourself.
The total effect :
The total effect is simply the sum of the two other effects so :
$$\frac {d\psi}{dt} = \frac {d\psi_l}{dt} + \frac {d\psi_s}{dt}$$
If you had multiple moons you would need multiple lunar correction terms.
Just for clarification that $\frac {d\psi}{dt}$ means the rate of change of the angle $\psi$ with respect to time $t$. To get how big an angle you'd move through in a century you do this :
$$\Delta \psi \approx \frac {d\psi}{dt} \Delta t$$
