I am not sure you need numerical simulation in this case. Aren't things solvable?
I guess one can think of the given "constraints" not so much as restrictions but as conservation laws. Basically, your equations of motion are the linear wave equations, so they have these "constraints" as first integrals. Therefore, if the initial conditions satisfy them at time $t=0$ then they continue to satisfies them for all $t$.
Basically, if you are looking for solutions of the form $\big(t,\,\vec{x}(t,\,\sigma)\big) \in \mathbb{R}^{1,D-1}$, where $(t,\sigma) \in \mathbb{R}^{1,1}$, change the independent time-space coordinates $(t,\sigma)$ into light coordinates:
\begin{align}
&t = \frac{u + v}{2} \\
&\sigma = \frac{u - v}{2}
\end{align}
Then your function $\vec{x} = \vec{x}\left(\frac{u + v}{2}, \, \frac{u - v}{2}\right)$ and the full space-time vector is
$$\left(\, \frac{u + v}{2}, \, \vec{x}\left(\frac{u + v}{2}, \, \frac{u - v}{2}\right)\,\right)$$
Observe that by the chain rule $${\partial}_v \,\vec{x} = {\partial}_t \,\vec{x} - {\partial}_{\sigma} \, \vec{x}$$ and hence
$${\partial}^2_{uv} \,\vec{x} = {\partial}^2_{tt} \,\vec{x} + {\partial}^2_{\sigma t} \,\vec{x} - {\partial}_{t \sigma}^2 \,\vec{x} - {\partial}_{\sigma \sigma}^2 \,\vec{x} = {\partial}^2_{tt} \,\vec{x} - {\partial}_{\sigma \sigma}^2 \,\vec{x} = 0$$
Thus, a vector function with the property that $\partial_{uv}\, \vec{x} = 0$ must have the form
$$\vec{x}(u,\,v) = \vec{x}_1(u) + \vec{x}_2(v)$$ for any single variable vector functions $\vec{x}_1(s)$ and $\vec{x}_2(s)$, with $s \in \mathbb{R}$, which means that when we switch back to the original variables
$$\vec{x}(t,\sigma) = \vec{x}_1(t + \sigma) \, + \, \vec{x}_2(t - \sigma)$$ Now we can check what the constraints impose:
\begin{align}
0 = \partial_t\, \vec{x} \cdot \partial_{\sigma} \, \vec{x} &= \big(\partial_t \, \vec{x}_1 + \partial_t\, \vec{x}_2\big) \cdot \big(\partial_{\sigma}\, \vec{x}_1 - \partial_{\sigma}\, \vec{x}_2\big) \\
&= \left(\frac{d\vec{x}_1}{ds}(t + \sigma) + \frac{d\vec{x}_2}{ds}(t - \sigma)\right) \cdot \left(\frac{d\vec{x}_1}{ds}(t + \sigma) - \frac{d\vec{x}_2}{ds}(t - \sigma)\right) \\
&= \left(\frac{d\vec{x}_1}{ds}(t + \sigma)\right)^2 - \left(\frac{d\vec{x}_2}{ds}(t - \sigma)\right)^2
\end{align} The latter should be true for any $(t,\sigma)$, so it should be true in particular for $(t = s, \sigma = 0)$, which is equivalent to
$$\left(\frac{d\vec{x}_1}{ds} (s) \right)^2 = \left(\frac{d\vec{x}_2}{ds}(s) \right)^2$$ and it should be true for $(t = 0, \sigma = s)$:
$$\left(\frac{d\vec{x}_1}{ds} (s) \right)^2 = \left(\frac{d\vec{x}_2}{ds}( - s)\right)^2 = \left(\frac{d\vec{x}_2}{ds}(s)\right)^2 = \left(\frac{d\vec{x}_1}{ds}( - s)\right)^2 $$
If you argue in an analogous way about the second constraint, you get
\begin{align}
R^2 = (\partial_t\, \vec{x})^2 + (\partial_{\sigma} \, \vec{x})^2 &= \big(\partial_t \, \vec{x}_1 + \partial_t\, \vec{x}_2\big)^2 + \big(\partial_{\sigma}\, \vec{x}_1 - \partial_{\sigma}\, \vec{x}_2\big)^2 \\
&= \left(\frac{d\vec{x}_1}{ds}(t + \sigma) + \frac{d\vec{x}_2}{ds}(t - \sigma)\right)^2 + \left(\frac{d\vec{x}_1}{ds}(t + \sigma) - \frac{d\vec{x}_2}{ds}(t - \sigma)\right)^2 \\
&= 2 \left(\frac{d\vec{x}_1}{ds}(t + \sigma)\right)^2 + 2\left(\frac{d\vec{x}_2}{ds}(t - \sigma)\right)^2
\end{align}
As before, the latter should be true for any $(t,\sigma)$, so it should be true in particular for $(t = s, \sigma = 0)$ and for $(t = 0, \sigma = s)$, which is equivalent to the conditions
$$\left(\frac{d\vec{x}_1}{ds}(s)\right)^2 + \left(\frac{d\vec{x}_2}{ds}(s)\right)^2 = \frac{R}{2}^2$$
$$\left(\frac{d\vec{x}_1}{ds}(s)\right)^2 + \left(\frac{d\vec{x}_2}{ds}(-s)\right)^2 = \frac{R}{2}^2$$ When you combine all of these conditions, you obtain the constraint
$$\left(\frac{d\vec{x}_1}{ds}(s)\right)^2 = \left(\frac{d\vec{x}_2}{ds}(s)\right)^2 = \frac{R^2}{4}$$
So if you pick any pair of vector functions $\vec{x}_1 : \mathbb{R} \to \mathbb{R}^{D-1}$ and $\vec{x}_2 : \mathbb{R} \to \mathbb{R}^{D-1}$ such that $$\left(\frac{d\vec{x}_1}{ds}(s)\right)^2 = \left(\frac{d\vec{x}_2}{ds}(s)\right)^2 = \frac{R^2}{4}$$ you can cosntruct the solution $$\vec{x}(t, \sigma) = \vec{x}_1(t + \sigma) + \vec{x}_2(t - \sigma)$$ which will satisfy the two constraints for any $(t, \sigma)$. For example, if you take any pair of vector functions $\vec{x}_j(s), \,\,\, j = 1,2$, describing a pair of curves in space $\mathbb{R}^{D-1}$ , parametrized by unit arc-length with respect to the Minkowski metric, restricted to the space-slices, then
$$\vec{x}(t, \sigma) = \vec{x}_1\left(\frac{R}{2}(t + \sigma)\right) + \vec{x}_2\left(\frac{R}{2}(t - \sigma)\right)$$ is a solution.