\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 240, pp. 1--12.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/240\hfil quasiperiodic self-oscillations] {Sensitivity analysis of stochastically forced quasiperiodic self-oscillations} \author[I. Bashkirtseva, L. Ryashko \hfil EJDE-2016/240\hfilneg] {Irina Bashkirtseva, Lev Ryashko} \address{Irina Bashkirtseva \newline Ural Federal University, Pr. Lenina 51, Ekaterinburg 620083, Russia} \email{irina.bashkirtseva@urfu.ru} \address{Lev Ryashko \newline Ural Federal University, Pr. Lenina 51, Ekaterinburg 620083, Russia} \email{lev.ryashko@urfu.ru} \thanks{Submitted May 30, 2016. Published August 31, 2016.} \subjclass[2010]{34D35, 37H10, 34K14} \keywords{Stochastic differential equations; tori; stochastic sensitivity; \hfill\break\indent quasipotential} \begin{abstract} We study a problem of stochastically forced quasi-periodic self-oscillations of nonlinear dynamic systems, which are modelled by an invariant torus in the phase space. For weak noise, an asymptotic of the stationary distribution of random trajectories is studied using the quasipotential. For the constructive analysis of a probabilistic distribution near a torus, we use a quadratic approximation of the quasipotential. A parametric description of this approximation is based on the stochastic sensitivity functions (SSF) technique. Using this technique, we create a new mathematical method for the probabilistic analysis of stochastic flows near the torus. The construction of SSF is reduced to a boundary value problem for a linear differential matrix equation. For the case of the two-torus in the three-dimensional space, a constructive solution of this problem is given. Our theoretical results are illustrated with an example. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{definition}[theorem]{Definition} \newtheorem{proposition}[theorem]{Proposition} \allowdisplaybreaks \section{Introduction} It is well known that complex multifrequency self-oscillations are basic operating modes for many real systems. Recent studies have shown that the variety observed in the behavior of nonlinear dynamical models under transition from order to chaos, is connected with a chain of bifurcations: equilibrium - periodic oscillations - quasiperiodic regime - chaos. Mathematically, each such transition is explained by the loss of stability of a simple attractor and the birth of a new, more complicated stable attractor. Here, an invariant manifold is a general mathematical model. Limit cycle is a corresponding invariant manifold for periodic oscillations. Quasi-periodic oscillations are modeled by invariant tori. Stability analysis of invariant manifolds and the study of their response on various perturbations (stochastic, for instance) is a key point in the understanding of the basic mechanisms of complex phenomena in real dynamical processes. An interplay of nonlinearity and stochasticity attract attention of many researchers \cite{Hor,Tel,LG}. An influence of random disturbances is widely studied both for ordinary differential equations \cite{ArnL,Has,Klod}, and for evolution and fractional dynamic systems \cite{Bez,Pra,Fla,Gess}. While stochastic cycles and equilibria are well-investigated objects of modern stability theory, appropriate analysis of randomly forced toroidal invariant manifolds is far from complete. After the investigations of Poincare, Denjoy and Arnold, the toroidal invariant manifolds are classical object of the modern qualitative theory of the deterministic differential equations \cite{Arn,Gu,Hol,Kat,Sam}. New nonlinear phenomena connected with torus canards are actively investigated in \cite{Ben,Kramer}. A probabilistic analysis of the noise-induced regimes in nonlinear systems with complex mixed-mode oscillations requires a study of corresponding stochastic attractors. A comprehensive mathematical description of stochastic attractors is given by the Kolmogorov-Fokker-Planck equation \cite{Gard}. However, for systems with dimensions two and higher, a direct use of these equations is very difficult. To avoid this complexity, various asymptotic and numerical approximations can be considered \cite{KS,Str}. For a weak noise, asymptotics based on the quasipotential function are widely used \cite{FW}. The quasipotential is a solution of the corresponding Hamilton-Jacobi equation. A development of the quasipotential approach has resulted in the elaboration of the theory for viscosity solutions of the Hamilton-Jacobi equation \cite{Cra,Lio,Per,Sch}. The quasipotential method continues to be actively developed (see, e.g. \cite{Far,Ish,Lo}). In many cases, it is important to describe a behavior of the stochastic system near the attractor of the initial deterministic system. A constructive method for the approximation of random states dispersion near limit cycles of nonlinear stochastic systems has been proposed in \cite{BR2004,RBGS09}. This method is based on the quadratic approximation of the quasipotential with the help of stochastic sensitivity function (SSF). SSF technique has been successfully used for the analysis of noise-induced transitions between attractors (limit cycles and equilibria) in \cite{BR11}, and for the suppression of noise-induced chaos in \cite{BCR12}. An extension of the SSF theory to the stochastic systems with toroidal attractors is an open novel problem even in the three-dimensional case. The aim of our work is to elaborate the stochastic sensitivity functions method for the analysis of randomly forced two-torus. In Section 2, for the constructive analysis of probabilistic distribution of randomly forced quasiperiodic oscillations, we present a new approach using a quadratic approximation of the quasipotential. A parametric description of this approximation is based on the SSF technique. Using this technique, we create a new mathematical model which specifies a dispersion of stochastic flows near the torus. The corresponding stochastic sensitivity matrix is a solution of the linear differential matrix equation. For the two-torus of a nonlinear system of stochastic differential equations, a parametric description of stochastic sensitivity function is derived. In Section 3, a case of the two-torus in three-dimensional space is studied in details. Here, SSF is scalar, and its construction is reduced to the solution of some functional equation. For the calculation of the stochastic sensitivity of two-torus, a constructive computer-oriented algorithm is suggested. The general theoretical results are illustrated on the example in Section 4, wherein a detailed parametric analysis of stochastic sensitivity of two-torus is given. \section{Stochastic sensitivity of torus} Consider the system \begin{equation} \dot{x}=f(x), \label{1} \end{equation} where $ x$ is $n $-vector, $ f $ is a sufficiently smooth vector-function. It is supposed that \eqref{1} has an invariant two-dimensional toroidal manifold $\Gamma$. Our analysis is based on the following parametrization of 2-torus $\Gamma$. Suppose that some closed sufficiently smooth curve $\alpha$ (equator) lies on the $\Gamma$ (see Figure \ref{fig1}). \begin{figure}[ht] \begin{center} \includegraphics[width=0.4\textwidth]{fig1} % ris1.eps \end{center} \caption{The line $\alpha$ is a closed curve (equator), $a=x(0,s)=\alpha(s)$ is an initial point of the solution $x(t,s)$, $b=x(T(s),s)=\alpha(\tau(s))$ is the first return point of the solution $x(t,s)$ on the curve $\alpha$.} \label{fig1} \end{figure} This curve is defined by the function $\alpha(s)$ on the interval $0 \leq s\leq 1$ with the condition $\alpha(0)=\alpha(1)$. Consider a solution $x(t,s)$ of \eqref{1} with the initial condition $x(0,s)=\alpha(s)$. It is supposed that the solution $x(t,s)$ leaves the point $\alpha(s)$ of the curve $\alpha$ and after rotation around the torus crosses the curve $\alpha$ again. Let $T(s)=\min\{t>0\;|\;x(t,s)\in\alpha\}$ be a first return time of the trajectory $x(t,s)$ on the curve $\alpha$ and $x(T(s),s)$ be a first return point. Let $\tau(s)$ be a point of the interval $[0,1)$ where $\alpha(\tau(s))=x(T(s),s)$. Here $\tau(s)$ is a Poincare first return function for intersections of the curve $\alpha$ by the phase trajectories of system \eqref{1}. A natural range of the variable $s $ is a circle. Consideration of the function $ \tau (s) $ on the interval $[0,1)$ only leads to discontinuity. To provide a continuity of the function $ \tau (s) $ we extend a domain of its definition on the infinite interval $- \infty < s < + \infty $. The equalities $x(t,s+1)=x(t,s)$, $x(T(s)+t,s)=x(t,\tau(s))$ allow us to extend a function $x (t, s)$ to the whole plane $ \Pi = \{(t, s) : -\infty < t < + \infty, \; -\infty < s < + \infty \}$. So, the function $x(t,s)$ defines one-to-one correspondence between points of 2-torus $\Gamma$ and points of the set $D=\{(t,s):0 \leq t < T(s),\;0\leq s<1\}$. The vector-functions $ \frac{\partial x(t,s)}{\partial t}$, $\frac{\partial x(t,s)}{\partial s} $ are linearly independent. For any point $\gamma\in\Gamma$ one can find $t=t(\gamma)$, $s=s(\gamma)$ such that $x(t,s)=\gamma$. In a neighbourhood $D$ of the torus $\Gamma$, consider a function $\gamma(x)$, where $\gamma(x)$ is a nearest point of 2-torus $\Gamma$ to $x$, ${ \Delta}(x)=x-\gamma(x)$ is a vector of a deviation of the point $x$ from the torus $\Gamma$. It is assumed that a neighbourhood $D$ is invariant for system \eqref{1}. \begin{definition} \label{def1} \rm A torus $\Gamma$ is called exponentially stable ($E$-stable) for system \eqref{1} in $D$ if there exist constants $K>0$, $l>0$ such that $$ \|{ \Delta}(x(t))\|\leq K e^{-lt} \|{ \Delta}(x_0)\|, $$ where $x(t)$ is a solution of \eqref{1} with initial condition $x(0)=x_0 \in D$. \end{definition} Stability of toroidal manifolds for deterministic systems was studied in \cite{Kol,Sam}. We consider the stochastic system \begin{equation} \dot x=f(x)+\varepsilon\sigma(x)\dot w, \label{2} \end{equation} where $w(t)$ is an $n$-dimensional Wiener process, $\sigma(x)$ is sufficiently smooth $n \times n$-function, $\varepsilon$ is a noise intensity. It is supposed that noise is nonsingular on the torus $\Gamma$ ($\det \sigma(x)|_\Gamma\neq 0$). Under stochastic disturbances random trajectories of \eqref{2} leave the torus $\Gamma$ and form some probability distribution around it. A general theoretical approach based on the Fokker-Planck-Kolmogorov (FPK) equation gives a most detailed description of such probabilistic distribution. However, a direct use of this equation is technically difficult, even in two-dimensional case. Therefore, asymptotic methods and approximations are commonly used \cite{Str,FW}. For weak noise, an asymptotic of the stationary probability density $\rho(x,\varepsilon)$ can be written in the form $$ \rho(x,\varepsilon)\approx K e^{-\frac{v(x)}{\varepsilon^2}}, $$ where the function $v(x)$ is a quasipotential \cite{FW}. The quasipotential being connected with some variational problem of the minimization of action functional is governed by the appropriate Hamilton-Jacobi equation. However, to solve this equation analytically is difficult too. In what follows, we will use a local description of the quasipotential in a small neighborhood of the torus $\Gamma$. The following function $v(x)$ is a quasipotential of \eqref{2} in the neighborhood $D$ of the torus $\Gamma$: $$ v(x)=\inf_{u \in U} J(x,u),\quad J(x,u)=\frac{1}{2}\int_0^\infty u^\top(\tau)\mathcal{S}^{-1}(y(\tau))u(\tau)d\tau,\; \mathcal{S}(y)=\sigma(y)\sigma^\top(y). $$ Here, $J(x,u)$ is an action functional, and $U$ is a set of admissible $n$-dimensional control inputs $u(\tau),\;0\leq\tau<+\infty$. These inputs $u\in U$ transfer a solution of the system \begin{equation} \dot y=-f(y)+u \label{3} \end{equation} from the initial point $y(0)=x$ to the torus $\Gamma : \lim_{\tau\to +\infty} { \Delta}(y(\tau))=0$. The function $v(x)$ is a solution of Hamilton-Jacobi equation \begin{equation} \big(f(x),\frac{\partial v}{\partial x}\big) + \frac{1}{2}\big(\frac{\partial v}{\partial x},\mathcal{S}(x)\frac{\partial v}{\partial x}\big)=0 \label{4} \end{equation} with the conditions \begin{equation} v|_\Gamma = 0, \quad ,v|_{D\setminus\Gamma} > 0. \label{5} \end{equation} From these properties, the quadratic form $$ \varphi (x) = \frac{1}{2}({ \Delta} (x) , \Psi (\gamma (x)){ \Delta} (x)) $$ is the first approximation of the quasipotential $v(x)$ near the torus $\Gamma: v(x)=\varphi(x)+O(\|{ \Delta}(x)\|^3)$. This quadratic form is parametrized by the symmetric non-negative $n\times n$-matrix function $\Psi(\gamma)= \frac{\partial^2 v}{\partial x^2}(\gamma)$ defined on the torus $\Gamma$. The matrix function $\Psi(\gamma)$ is singular. Indeed, consider an arbitrary smooth curve $\beta$ lying on $\Gamma$ and passing through a point $\gamma\in\Gamma$. Let this curve have the parametrization: $\beta(\tau)$, $\beta(0)=\gamma$. Differentiating the identity $$\frac{\partial v}{\partial x}(\beta(\tau))\equiv 0$$ with respect to $\tau$, we have $$ \frac{\partial^2 v}{\partial x^2}(\beta(\tau))\frac{d\beta(\tau)}{d\tau} \equiv 0. $$ This means that \begin{equation} \Psi(\gamma)r(\gamma)\equiv 0, \label{6} \end{equation} where $r(\gamma)$ is an arbitrary vector tangent to $\Gamma$ at the point $\gamma$. The condition \eqref{6} means that the matrix $\Psi(\gamma)$ is singular and it holds that $\operatorname{rank}(\Psi(\gamma))\leq n-2$. For the function $\Psi(\gamma)$, consider the following parametrization. Using a family of solutions $x(t,s)\in\Gamma$ of the deterministic system \eqref{1}, one can introduce the matrix function \begin{equation} V(t,s)=\Psi(x(t,s)) \label{7} \end{equation} defined on $\Pi$. It follows from $$ x(t,s+1)=x(t,s),\;x(T(s)+t,s)=x(t,\tau(s)), $$ that \begin{equation} V(t,s+1)=V(t,s), \quad V(T(s)+t,s)=V(t,\tau(s)). \label{8} \end{equation} Note that the singularity condition \eqref{6} implies $$ V(t,s)\frac{\partial x(t,s)}{\partial t} \equiv 0,\quad V(t,s)\frac{\partial x(t,s)}{\partial s} \equiv 0. $$ This system can be rewritten in the form \begin{equation} V(t,s)y(t,s)=0, \quad V(t,s)z(t,s)=0. \label{9} \end{equation} Here, \begin{equation} y(t,s)=f(x(t,s)), \label{10} \end{equation} and $z(t,s)$ is a solution of the system \begin{equation} \frac{\partial z}{\partial t}=F(t,s)z\,,\quad F(t,s)=\frac{\partial f}{\partial x}(x(t,s)) \label{11} \end{equation} with the initial value $z(0,s)= \frac{d\alpha(s)}{ds}$, Consider a space $\Sigma$ of the sufficiently smooth symmetric $n\times n$-matrix functions $V(t,s)$ defined on the plane $\Pi$ with conditions \eqref{8}, \eqref{9}. Let $P_{y,z}$ be a matrix of the projection onto subspace that is orthogonal to the vectors $y$ and $z$. Consider the matrix function $P(t,s)=P_{y(t,s),z(t,s)}$, where $y(t,s)$ and $z(t,s)$ are defined in \eqref{10}, \eqref{11}. Note that $$ P_{y,z}=P_y - \frac{P_yzz^{\top} P_y}{z^{\top} P_yz} \,, \quad P_y=I-\frac{yy^{\top} }{y^{\top} y} \,,\quad \operatorname{rank}(P_{y,z})=n-2. $$ \begin{definition} \label{def2} \rm A matrix $V(t,s)\in\Sigma$ is called by $P(t,s)$-positive definite at the point $(t,s)$ if for any vector $u$ satisfying $P(t,s)u\neq 0$ the inequality $(u,V(t,s)u)>0$ holds. The matrix $V(t,s)$ which is $P(t,s)$-positive definite for any $(t,s)\in\Pi$, is called by $P$-positive definite. \end{definition} In the space $\Sigma$, consider a cone $\mathcal{K}=\{V\in\Sigma|V(t,s)$ is non-negative definite for any $(t,s)\in\Pi \}$ and a set $\mathcal{K}_P=\{V\in\Sigma : V\text{ is a $P$-positive definite matrix}\}$. Taking into account a sufficient smoothness of the quasipotential $v(x)$, one can differentiate the Hamilton-Jacobi equation \eqref{4} and substitute $x=x(t,s)$. As a result, for $ V(t,s)$, we obtain the matrix differential Bernoulli equation \begin{equation} \frac{\partial V(t,s)}{\partial t}+F^\top(t,s)V(t,s)+V(t,s)F(t,s) +V(t,s) S(t,s) V(t,s) =0, \label{12} \end{equation} for $V\in\Sigma$, where $$ F(t,s)=\frac{\partial f(x(t,s))}{\partial x},\quad S(t,s) = \mathcal{S}(x(t,s)) =\sigma (x(t,s))\sigma^\top(x(t,s)). $$ Along with \eqref{12} consider the linear matrix equation \begin{equation} \frac{\partial W(t,s)}{\partial t}=F(t,s)W(t,s)+ W(t,s)F^\top(t,s) + P(t,s)S(t,s)P(t,s), \label{13} \end{equation} for $W\in\Sigma$. \begin{proposition} \label{prop1} Let the torus $\Gamma$ of system \eqref{1} be $E$-stable. Then \eqref{13} has a unique solution $W\in\mathcal{K}$. \end{proposition} \begin{proof} The $E$-stability of the torus means \cite{Non2004} that for any element $C\in\mathcal{K}$ there exists a unique solution $V\in\mathcal{K}$ of the matrix Lyapunov equation $L[V]=-C$ with the following Lyapunov operator $$ L[V]= \frac{\partial V(t,s)}{\partial t}+F^\top(t,s)V(t,s)+V(t,s)F(t,s). $$ Since the cone $\mathcal{K}$ is reproducing in $\Sigma$ \cite{Krasn}, the operator $L$ is invertible on the entire space $\Sigma$. For the space $\Sigma$ with the scalar product $$ \langle V,W\rangle=\lim_{T \to \infty}\frac{1}{T} \int_{0}^{T} \operatorname{tr}(V(t,s)W(t,s))dt, $$ the conjugate operator of $L$, has the representation $$ L^*[W]= -\frac{\partial W(t,s)}{\partial t}+F(t,s)W(t,s)+W(t,s)F^\top(t,s). $$ Then \eqref{13} can be written as $L^*[W]+PSP=0$. The matrix $ W=-(L^{-1})^*[PSP] \in \mathcal{K}$ is a unique solution of this equation. The proof is complete. \end{proof} \begin{proposition} \label{prop2} Let $W$ be a solution of \eqref{13}. Then, the the matrix $V=W^+$ is a solution of \eqref{12}. \end{proposition} \begin{proof} Multiplying \eqref{13} by $V$ from the left and right, we obtain \begin{equation} V\dot W V=VFP+PF^\top V+VSV. \label{14} \end{equation} Here, $\dot W = \frac{\partial W}{\partial t}$. Differentiating the equality $V=VWV$, we obtain $\dot V=\dot V P+V \dot W V+P \dot V$, and $V \dot W V =\dot V-\dot V P- P\dot V$. Multiplying the last equality by $P$ from the left and right, one can obtain $V \dot W V=-P\dot V P$. Taking into account \eqref{14} we have \begin{equation} P [\dot V+VF+F^\top V+VSV]P=0. \label{15} \end{equation} Consider a vector $r=r(t,s)$ that belongs to the plain spanned by the vectors $y(t,s), z(t,s)$. Since $V r=0$, we obtain \begin{equation} [\dot V+VF+F^\top V+VSV]r=[\dot V+VF] r =\dot{[Vr]}=0. \label{16} \end{equation} It follows from \eqref{15}, \eqref{16} that $V=W^+$ is a solution of \eqref{12}. The proof is complete. \end{proof} Using these results, one can write the following approximation of the quasipotential: \begin{equation} v(x) = \varphi(x) +O(\|{ \Delta} (x)\|^3),\label{17} \end{equation} where $$ \varphi(x)= \frac{1}{2} ({ \Delta} (x) , \Phi^+ (\gamma (x)){ \Delta} (x)),\;\;\Phi (\gamma)=W(t(\gamma),s(\gamma)). $$ Here, $W(t,s)$ is a solution of the matrix equation \eqref{13}. This approximation $\varphi(x)$ of the quasipotential $v(x)$ allows us to represent an asymptotic of the stationary density in the form of the Gaussian distribution: $$ \rho(x,\varepsilon)\approx \rho_\ast (x,\varepsilon) = K \exp\Big(- \frac{({ \Delta} (x) , \Phi^+ (\gamma (x)){ \Delta} (x))}{2 \varepsilon^2}\Big) $$ with covariance matrix $\varepsilon^2\Phi(\gamma)$. This matrix characterizes a dispersion of random trajectories of \eqref{2} at the point $\gamma$ of the toroidal surface. Let $\lambda_1(\gamma)\geq\lambda_2(\gamma)\geq \dots \geq\lambda_n(\gamma)$ be eigenvalues, and $h_1(\gamma)$, $h_2(\gamma),\dots ,h_n(\gamma)$ be orthonormal basis of eigenvectors of the matrix $\Phi(\gamma)$. Since for any $\gamma\in\Gamma$ the matrix $\Phi(\gamma)$ is singular ($\operatorname{rank}\Phi(\gamma)=n-2$) then $\lambda_{n-1}(\gamma)=\lambda_n(\gamma)\equiv 0$. Corresponding eigenvectors $h_{n-1}(\gamma)$ and $h_n(\gamma)$ belong to the plane that is tangent to torus. For noises non-degenerate on $\Gamma$, other eigenvalues $\lambda_1,\dots ,\lambda_{n-2}$ are positive and define a dispersion of random trajectories in direction of the vectors $h_1, \dots , h_{n-2}$. These vectors specify a basis of the normal hyperplane. The matrix $\Phi(\gamma)$ characterizes a response of system \eqref{2} near the torus $\Gamma$ on the random input. If we consider \eqref{2} as a converter of the stochastic input (stationary Wiener process $w(t)$) to the stochastic output (stationary distribution of random trajectories near $\Gamma$), the eigenvalues of $\Phi(\gamma)$ define coefficients of amplification $(\lambda_i>1)$ and weakening $(\lambda_i<1)$ of this converter. The stochastic sensitivity function $\Phi(\gamma)$ allows us to describe the nonuniformity of the dispersion of random trajectories near torus along all directions, and to mark more and less sensitive to noise portions of the toroidal surface. \section{Stochastic sensitivity of a 2-torus in the three-dimensional space} Consider system \eqref{2} for $n=3$ In this case, the projective matrix $P(t,s)$ has a rank equal to one, and so this matrix can be written as $P(t,s)=p(t,s)p^\top(t,s)$ where $p(t,s)$ is a unit vector orthogonal to vectors $y(t,s)$ and $z(t,s)$. The matrix $W$ which is a solution of \eqref{13} characterizes a dispersion of random trajectories near 2-torus. This matrix also has a rank equal to one, and can be written as $W(t,s)=\mu(t,s)P(t,s)$. The projective matrix $P(t,s)$ can be found uniquely for any point of the initial deterministic surface $\Gamma$. So, a problem of the construction of the stochastic sensitivity function is reduced here to the calculation of the scalar function $\mu(t,s)$. This function describes a dispersion of random trajectories in the direction that is orthogonal to the torus $\Gamma$ at the point $x(t,s)$. \begin{proposition} \label{prop3} The matrix $W(t,s)=\mu(t,s)P(t,s)$ is a solution of \eqref{13} if and only if the scalar function $\mu(t,s)$ is a solution of the system \begin{equation} \frac{\partial \mu}{\partial t}(t,s) =a(t,s)\mu(t,s)+b(t,s) \label{18} \end{equation} with the conditions \begin{gather} \mu(t,s+1)=\mu(t,s), \label{19} \\ \mu(T(s)+t,s)=\mu(t,\tau(s)). \label{20} \end{gather} Here, $$ a(t,s)=p^\top(t,s)[F^\top(t,s)+F(t,s)]p(t,s), \quad b(t,s)=p^\top(t,s)S(t,s)p(t,s). $$ \end{proposition} \begin{proof} Substituting the solution $W(t,s)=\mu(t,s)p(t,s)p^\top(t,s)$ into the equation \eqref{13}, we obtain $$ \dot\mu p p^\top + \mu (\dot p p^\top+p (\dot p)^\top)= \mu(F p p^\top +p p^\top F^\top)+p p^\top S p p^\top. $$ Multiplying this equality by $ p^\top$ from the left and by $p$ from the right, and taking into account that $p^\top p\equiv 1$, $\dot{(p^\top p)}=\dot{p}^\top p+p^\top \dot{p}\equiv 0$, we obtain the following equation for $\mu$: $$ \dot{\mu}=p^\top (F+F^\top) p \mu +p^\top S p. $$ A proof of the converse implication follows from calculations written above. The proof is complete. \end{proof} As a result, for $n=3$, the stochastic sensitivity function is found from the scalar equation \eqref{18} with conditions \eqref{19}, \eqref{20}. A general solution of \eqref{18} is as follows \begin{equation} \mu(t,s)=g(t,s)[c(s)+h(t,s)] \label{21} \end{equation} where functions $g(t,s)$ and $h(t,s)$ can be written in the explicit form \begin{equation} g(t,s)=\exp\Big(\int_0^t a(\tau,s)d\tau\Big),\quad h(t,s)=\int_0^t \frac{b(\tau,s)}{g(\tau,s)}d\tau. \label{22} \end{equation} The unknown function $c(s)$ plays a role of initial data for $\mu(t,s)$: $\mu(0,s)=c(s)$. Condition \eqref{19} gives the equation $c(s+1)=c(s)$. Condition \eqref{20} implies $$ g(T(s),s)[c(s)+h(T(s),s)]=c(\tau(s)). $$ Denote \begin{equation} \alpha(s)=g(T(s),s),\quad \beta(s)=\alpha(s)h(T(s),s). \label{23} \end{equation} As a result, for the function $c(s)$ we obtain a functional equation \begin{equation} c(\tau(s))=\alpha(s)c(s)+\beta(s) \label{24} \end{equation} with the condition \begin{equation}c(s+1)=c(s). \label{25} \end{equation} So, a construction of the stochastic sensitivity function $\mu(t,s)$ for 2-torus in three-dimensional space is reduced to the solution of the functional equation \eqref{24} with condition \eqref{25}. For $E$-stable torus $\Gamma$, an existence and uniqueness of the solution of \eqref{24}, \eqref{25} directly follows from the Propositions \ref{prop1} and \ref{prop3}. A solution $\bar c(s)$ of \eqref{24},\eqref{25} can be found by the stabilization method. Consider sequences $ s_0=s, s_1 ,\dots , s_k,\dots $, where $ s_{k+1}=\tau (s_k)$, and $\bar c_0,\bar c_1 ,\dots ,\bar c_k,\dots $, where $ \bar c_k=\bar c(s_k)$. Because of \eqref{24}, the values $ \bar c_k$ are connected by the equation \begin{equation} \bar c_{k+1} = \alpha_k \bar c_k +\beta_k , \quad \alpha_k= \alpha (s_k), \quad \beta_k =\beta (s_k). \label{26} \end{equation} For elements $ \bar c_k $, consider approximations $ c_k$ generated by the recurrence formula $$ c_{k+1} = \alpha_k c_k +\beta_k, $$ where $ c_0 $ is some approximation for $ \bar c_0$. \begin{proposition} \label{prop4} Let the torus $\Gamma$ be $E$-stable. Then $\lim_{k \to \infty}( c_k - \bar c_k)=0$ for any initial approximation $ c_0$. \end{proposition} \begin{proof} The error $r_k= c_k - \bar c_k$ satisfies the equation $r_{k+1} = \alpha_k r_k $, and relates with the initial error $ r_0 $ by the formula $r_k = q_k r_0$, where $$ q_k=\prod_{i=0}^{k-1}\alpha_i=\exp\Big(\int_0^{T_k(s)}a(t,s)dt\Big). $$ Here, $T_k(s)$ is passing time of the trajectory $x(t,s)$ along the scroll consisting of $k$ turns. \end{proof} A necessary and sufficient condition of the $E$-stability of the torus in system \eqref{1} is the following \cite{Non2004}: \begin{equation} \max_s\lim_{T\to\infty}\frac{1}{T}\int_0^{T}a(t,s)dt<0. \label{27} \end{equation} Condition \eqref{27} implies that $\lim_{k\to\infty} q_k=0$, and therefore $\lim_{k\to\infty} r_k=0$ independently on the choice of $c_0$. These mathematical results can be summarized in the following computer-oriented algorithm. \subsection*{Algorithm for the calculation of stochastic sensitivity of 2-torus} \quad \begin{itemize} \item[(1)] Let $s_0, s_1,\dots, s_N$ be a discretization of the interval $0\leq s \leq 1$, where $s_{k+1}=\tau(s_k)$. For the each point $s_k$, find a solution $x(t,s_k)$ of the equation \eqref{1} with initial condition $x(0,s_k)=\alpha(s_k)$ on the interval $0\leq t \leq T(s_k)$. \item[(2)] Using \eqref{22} with \begin{gather*} a(t,s)=p^\top(t,s)[F^\top(t,s)+F(t,s)]p(t,s), \\ b(t,s)=p^\top(t,s)S(t,s)p(t,s),\\ F(t,s)=\frac{\partial f}{\partial x}(x(t,s)),\quad S(t,s)=\sigma(x(t,s))\sigma^\top(x(t,s)), \end{gather*} calculate functions $g(t,s_k),\;h(t,s_k)$. \item[(3)] Find $\alpha(s_k)$ and $R\beta(s_k)$ from formulas \eqref{23}. \item[(4)] By the iterative method \eqref{26}, find a solution $c(s)$ of the functional equation \eqref{24} with condition \eqref{25} at the points $s_k$. \item[(5)] Finally, find the stochastic sensitivity functions $\mu(t,s_k)$ by the formula \eqref{21}. \end{itemize} \section{Example} In the three-dimensional phase space of variables $(x,y,z)$, we consider the 2-torus $\Gamma$ defined by $$ (\sqrt{x^2+y^2}-1)^2+z^2=r_0^2,\quad 0