辛扭转映射的 KAM 迭代算法

在一个动力系统的相空间中,拟周期轨道一般是指这样的不变子集: 其在拓扑上与一个 $n$ 维环面同胚, 在动力学上, 往往存在着该环面上的某个参数化, 使得系统的动力学限制在环面上表现为简单的平移(针对映射系统)或者是线性流(针对流系统). 一般而言, 环面可能在系统的扰动下得以保持, 但环面上的流在扰动下可能会展现出极其复杂的性态.

在上述意义下, 拟周期轨道的存在性目前似乎只有唯一的研究工具: KAM 理论. 拟周期轨道的存在性往往可以归结为某个算子方程的解, KAM 理论通过线性化这个算子, 再利用频率的代数性质以及函数空间的性质, 往往可以求得这个线性化的算子的唯一解. 而后再精细地控制该解的定义区域以及残差项的误差, 通过迭代的方式得到原来非线性算子方程的精确解. 这个步骤包含了算子的线性化过程, 且一次迭代的残差项也是以上一次残差项的二阶小量的形式, 因此上述方法也称作是改进的牛顿迭代法. 几乎所有的 KAM 定理均涉及到了上述过程.

值得注意的是, 上述改进的牛顿迭代法完全可以在数值上进行实现, 并且实现的算法具有非常好的速度以及准确度, 同时该算法也能初步判断某个频率下环面是否破坏: 因为若该频率的拟周期解不存在, 算法是发散的. 下面我们就来详细介绍这种算法, 并具体地讨论其如何在一般的通有计算机语言上实现.

问题描述

我们采用的模型是标准映射:

$$\begin{aligned} x'&=x+y+\frac{k}{2\pi}\sin(2\pi x),\\ y'&=y+\frac{k}{2\pi}\sin(2\pi x). \end{aligned}$$

最后我们会简要说明如何把这个方法应用于某些特殊的微分方程, 这些微分方程的某个截面所定义的映射的生成函数必须是显式给出的. 这个条件似乎过于苛刻, 但这种方法仍然不失为一种强大的工具.

标准映射是一个辛扭转映射. 具体来说, 定义生成函数:

$$h(x,x')=\frac{1}{2}(x'-x)^2-\frac{k}{(2\pi)^2}\cos(2\pi x).$$

那么标准映射与下列方程等价:

$$\frac{\partial h}{\partial x}=-y,\ \frac{\partial h}{\partial x'}=y'.$$

因此, 若一序列 $\{x_{i}\}_{i\in\mathbb{Z}}$ 满足

$$h_2(x_{i-1},x_{i})+h_1(x_i,x_{i+1})=0, \forall i\in\mathbb{Z},$$

那么由 $y_{i}=h_{2}(x_{i-1},x_{i})$ 定义的 $\{(x_i,y_i)\}_{i\in\mathbb{Z}}$ 确实是标准映射的一条轨道.

实际上, 如果我们取标准映射相空间中的参数化曲线: $\Gamma(x)=(u(x), -h_{1}(u(x-\omega),u(x)))$. 那么不难验证有: $T(\Gamma(x))=\Gamma(x+\omega)$. 这符合我们在前言中对拟周期解的描述.

对于某个固定的 $\omega$, 下面我们就来讨论如何求解算子方程:

$$E(u):=h_2(u(x-\omega),u(x))+h_1(u(x),u(x+\omega))=0$$

的解.

算子方程的线性化

为了记号的方便, 我们令 $u^{+}=u(x+\omega),\ u^{-}=u(x-\omega),\ h_{ij}=h_{ij}(u,u^+),\ h_{ij}^{-}=(u^-,u).$

先前我们提到过, 求解算子方程 $E(u)=0$ 用到了改进的牛顿迭代法. 所以下面我们进行通常的牛顿迭代步骤. 首先对 $u$ 施加一个小的周期扰动 $v$, 这样 $u+v$ 仍然满足原先 $u$ 满足的条件. 那么直接对 $v$ 展开得到:

$$E(u+v)=E(u)+(h_{11}+h_{22}^{-})v+h_{12}v^{+}+h_{21}^{-}v^{-}+Q,$$

其中 $Q$ 代表了高阶项. 为了记号上的方便, 我们可以定义下面的算子:

$$E'(u)(v):=(h_{11}+h_{22}^{-})v+h_{12}v^{+}+h_{21}^{-}v^{-},$$

这里算子 $E'(u)$ 表示 $E(u)$ 在 $u$ 处的导数.

因此, 如果按照牛顿迭代的方法, 那么首先我们应当去求解合适的 $v$ 使得

$$E'(u)(v)=-E(u).$$

直接在数值上求解 $v$, 再令新的 $u$ 为 $\bar{u}=u+v$, 继续重复上述步骤的方法称为拟周期轨道的直接法. 我们这里简要叙述一下直接法求解的步骤, 可以看到其最终是求一个大矩阵的逆.

直接法

首先, 我们先要取定初始的 $u$ 才能进行后续的计算. 一般地, 直接取 $u$ 为平移即可, 即 $u(\theta)=\theta +\omega$. 我们将所要求解的式子展开为:

$$(h_{11}+h_{22}^{-})v+h_{12}v^{+}+h_{21}^{-}v^{-}=-h_2(u^-,u)-h_1(u,u^+)$$

在数值上, 表示一个函数一般用等分布地 $N$ 个点的函数值去近似一个函数. 由于我们所要处理的函数都是定义在 $[0,1]$ 上, 故我们不妨取这 $N$ 个点为 $\{0,\frac{1}{N},\frac{2}{N},...,\frac{N-2}{N},\frac{N-1}{N}\}$.

在一般的某一次迭代中, 我们仅仅知道函数 $u$ 在这 $N$ 个点上的分布. 要进行下一次迭代, 我们需要知道 $v$ 在这 $N$ 个点上的分布才行. 而上述方程的未知量似乎太多了: $u^-,u^+,v^-,v,v^+$. 实际上, 利用离散 Fourier 变换, 假设未知量是函数 $v$ 在 $N$ 个点的函数值, 那么 $v^-$ 以及 $v^+$ 在这些点的函数值可以由 $v$ 的函数值线性表出. 同样地, 可以利用离散 Fourier 变换求得 $u^{-},u^+$ 的分布, 那么上述方程中未知的仅有 $v$ 了. 由于函数值有 $N$ 个, 可以建立 $N$ 个有 $N$ 个未知数的方程. 注意到, 由于 $v^-$ 以及 $v^+$ 的函数值可以由 $v$ 的函数值线性表出, 这些方程是线性方程! 所以最后求 $v$ 的函数值的问题变成求矩阵逆的问题.

同调方程

引入两个差分算子:

$$\nabla f(\theta)=f(\theta+\omega)-f(\theta),\ \nabla^{*}f(\theta)=f(\theta)-f(\theta-\omega).$$

这里我们直接给出线性化算子:

$$\nabla^{*}(h_{12}u_\theta u_\theta^+ \nabla w)=-u_\theta E(u),$$

其中 $v=u_\theta w$, 迭代格式为 $\bar{u}=u+v$. 具体推导细节可参见文献 [1]. 线性化算子的未知量是 $w$. 其余量如 $u,u_\theta$ 由上一次迭代的信息给出.

离散 Fourier 变换及其简单性质

在上述讨论中, 我们提到了离散 Fourier 变换的概念, 下面我们来详细介绍这个变换.

假设有 $N$ 个数据点: $\{x_0,x_1,...,x_{N-1}\}$. 离散 Fourier 变换将 $N$ 个点变为令一组 $N$ 个点 $\{y_0,y_1,...,y_{N-1}\}$, 其中

$$y_k=\frac{1}{N}\sum_{j=0}^{N-1}x_{j}e^{-2\pi i\frac{jk}{N}}.$$

下面我们尝试解释这个变换的涵义. 我们可以将 $x_k$ 理解为一个函数 $f$ 在 $[0,1]$ 上的等分布函数值, 这些点为 $\{0,\frac{1}{N},\frac{2}{N},...,\frac{N-2}{N},\frac{N-1}{N}\}$. 那么可以看出, 上述求和恰好是积分

$$c_{k}=\int_{0}^{1}f(x)e^{-2\pi i k x}dx$$

的黎曼和. 这样直观地, 我们似乎可以认为 $\{y_0,y_1,...,y_{N-1}\}$ 近似地表示了函数 $f$ 的 Fourier 系数 $\{c_0,c_1,...,c_{N-1}\}$. 然而并非如此, 原因是当 $y_{k}$ 的指标 $k$ 太大时, $e^{-2\pi i k x}$ 的频率是 $\frac{1}{k}$, 即其震荡是相当快的, 这样用如此少的点进行黎曼求和会带来相当大的误差.

解决这个问题的方法是应用数列 $y_{k}$ 的周期性. 首先注意到, 在 $y_{k}$ 的定义中实际上并不需要 $k$ 一定在 $0,...,N-1$ 个数之间. 这样我们可扩充 $y_{k}$ 到整个整数指标集: $\{y_{k}\}_{k\in\mathbb{Z}}$. 那么我们有:

$$y_{N+k}=\frac{1}{N}\sum_{j=0}^{N-1}x_{j}e^{-2\pi i\frac{j(N+k)}{N}}=\frac{1}{N}\sum_{j=0}^{N-1}x_{j}e^{-2\pi i\frac{jk}{N}}=y_k.$$

这意味这序列 $y_{k}$ 是 $N$ 周期的. 对于 $k\in\{-\frac{N}{2},...,-1\}$, 显然 $y_{k}=y_{k+N}$ 对于 $c_{k}$ 的逼近要好于 $c_{k+N}$, 这是由于 $k$ 的绝对值的指标要比 $k+N$ 的绝对值的指标要小, 因此周期函数 $e^{-2\pi i kx}$ 的震荡要弱一些, 黎曼求和的精度也要高一些.

所以, 通常地, 对于 $N$ 是偶数, 我们一般用

$$\{y_0,y_1,...,y_{N-1}\}\approx\{c_0,c_1,...,c_{\frac{N}{2}-1},c_{-\frac{N}{2}},...,c_{-1}\}.$$

对于 $N$ 是奇数, 我们一般用

$$\{y_0,y_1,...,y_{N-1}\}\approx\{c_0,c_1,...,c_{\frac{N-1}{2}},c_{-\frac{N-1}{2}},...,c_{-1}\}.$$

实际上, 有一些文章专门讨论了离散 Fourier 变换与真实函数的精度问题. 这里我们不作讨论.

类似地, 可以定义离散 Fourier 逆变换:

$$x_k=\sum_{j=0}^{N-1}y_{j}e^{2\pi i\frac{jk}{N}}.$$

另外, 注意到可以用矩阵乘法来定义上述离散 Fourier 变换. 假设矩阵:

$$\mathbf{M}=\left(\begin{array}{ccccc} 1 & 1 & 1 & \cdots & 1 \\ 1 & e^{-2 i \pi 1 \times 1 / N} & e^{-2 i \pi 2 \times 1 / N} & \cdots & e^{-2 i \pi 1 \times(N-1) / N} \\ 1 & e^{-2 i \pi 1 \times 2 \times 1 / N} & e^{-2 i \pi 2 \times 2 / N} & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & e^{e-2 i \pi(N-2) \times(N-1) / N} \\ 1 & e^{-2 i \pi(N-1) \times 1 / N} & \cdots & e^{e-2 i \pi(N-1) \times(N-2) / N} & e^{-2 i \pi(N-1) \times(N-1) / N} \end{array}\right)$$

那么我们有:

$$y=\frac{1}{N} \mathbf{M}x.$$

这里我们总结一下离散 Fourier 变换及其逆变换的作用. 知道一个函数在区间上的等分布值, 我们可以利用离散 Fourier 变换求得其 Fourier 系数. 同样, 知道一个函数的 Fourier 系数, 并将其排列成上述特定的形式, 利用离散 Fourier 逆变换可以反求等分布的函数值.

数值上有一种快速计算 Fourier 变换的方法, 称作是 FFT (Fast Fourier Transformation). 这种方法可以大幅缩减离散 Fourier 变换及其逆变换的计算量. 在许多程序语言中均内置了 FFT, 如 Mathematica 的 Fourier 与 InverseFourier 函数, Matlab 中的 fft, ifft, fft2 等, Julia 的 FFTW 程序包等. 值得注意的是, 这些函数定义的方式有所不同, 例如 Matlab 里默认的 fft 是不除以列表的长度的, Mathematica 使用我们所定义的离散 Fourier 变换以及逆变换, 需要设置 Fourier 变换的参数为 {-1,-1}.

具体算法的实现

直接法

我们再来回顾一下直接法需要做的事: 假设已知 $u$ 在区间上的均匀分布, 我们需要求 $v$ 在区间上的均匀分布, 而这两者由如下方程联系起来:

$$(h_{11}+h_{22}^{-})v+h_{12}v^{+}+h_{21}^{-}v^{-}=-h_2(u^-,u)-h_1(u,u^+)$$

算法实现:

• 假设函数 $v$ 的等分布由 $N$ 个未知量表出: $\{v_0,...,v_{N-1}\}$.

• 计算函数 $v^-$ 的分布, 通过以下方法.

  • 首先利用离散 Fourier 变换求得 $v$ 的 Fourier 系数: $v_f=\{v_{0f},...,v_{-1f}\}$.

  • 再将 $v_f$ 的每个元素 $v_{jf}$ 乘以 $e^{-2\pi ij \omega }$, 得到 $v_{pf}$.

  • 再对 $v_{pf}$ 做离散 Fourier 逆变换, 便得到了 $v^-$ 的分布.

• 同样方法计算 $v^{+}$ 的分布, 不过进行 Fourier 变换后要乘以系数 $e^{2\pi ij \omega }$.

• 同样方法计算 $u^{-}, u^{+}$ 的分布.

• 在 $[0,1]$ 上的每个分布点线性化方程都成立. 对上述方程的每个分布点建立方程, 得到 $N$ 个方程.

• 由于只有 $\{v_0,...,v_{N-1}\}$ 个未知量. 求解上述线性方程组得到 $v$.

• 令新的 $u$ 为 $u+v$. 循环实现一次.

• 重复上述循环.

同调方程法

同样, 我们可以回顾一下牛顿迭代格式由同调方程给出:

$$\nabla^{*}(h_{12}u_\theta u_\theta^+ \nabla w)=-u_\theta E(u).$$

上述方程可以分为两步求解:

$$\begin{cases} &\nabla^{*}\psi=g,\\ &\nabla w= p(\psi+\mu), \end{cases}$$

其中 $p=(h_{12}u_\theta u_\theta^+)^{-1}$, $g=-u_\theta E(u)$,

$$\mu=-(\int_{0}^{1}p(\theta)d\theta)^{-1}\int_{0}^{1}p\psi d\theta.$$

算法实现

• 首先按照直接法中所叙述的方法, 计算 $u^+, u^-$ 的分布值.

• 下面计算 $u_{\theta}$ 的分布值:

  • 利用离散 Fourier 变换求得 $u$ 的 Fourier 系数 $u_{fi}$, $i\in\{0,...,-1\}$.

  • 对每个 $u_{fj}$ 乘以 $2\pi j$ 得到 $u_{\theta}$ 的 Fourier 系数

  • 利用 Fourier 逆变换求得 $u_{\theta}$ 的分布值.

• 类似地可以计算 $u_{\theta}^{+}$ 的分布值, 但乘法系数变为 $2\pi j e^{2\pi i j}$.

• 计算函数 $g$ 的分布值

• 利用 $g$ 的函数分布值求得其 Fourier 系数 $g_{fj}$.

• $g_{fj}$ 乘以 $\frac{1}{1-e^{-2\pi \omega i j}}$ 得到函数 $\psi$ 的 Fourier 系数.

• 使用 Fourier 逆变换得到函数 $\psi$ 的分布值.

• 计算 $p,p\psi$ 的分布值, 并使用黎曼求和计算 $\mu$.

• 计算 $p(\psi+\mu)$ 的分布值, 并利用 Fourier 变换得到 Fourier 系数.

• $p(\psi+\mu)$ 的 Fourier 系数 $pw_{fj}$ 乘以 $\frac{1}{e^{2\pi \omega i j}-1}$, 得到 $w$ 的 Fourier 系数.

• 对 $w$ 的 Fourier 系数做 Fourier 逆变换得到 $w$ 的函数值分布.

• 令新的 $u$ 为 $u+u_{\theta}w$, 完成一次循环.

• 重复循环.

讨论: 虽然同调方程方法涉及到的理论知识过多, 且程序实现起来也繁琐一些, 但其准确度以及速度是直接法无法比拟的. 在下面的 Julia 代码中, 我们取 $N$ 为 2^12=4096, $\omega=\frac{\sqrt{5}-1}{2}$, $k=0.6$, 初始的 $u$ 即为平移, 这样大量的数据下, 整个牛顿迭代过程完成 100 次迭代仅不到 0.3 秒. 而同样的直接法, 由于涉及到符号的运算以及求导运算, 取 N=100, 迭代几次则需要花费十几秒.

程序实现

同调方程 Julia 代码:

using FFTW
const n = 2^12
const ω = 2 \ (sqrt(5) - 1)
const k = 0.6

using Symbolics

@variables k s0 s1

hh = 0.5 * (s1 - s0)^2 - ((2 * π)^2 \ k) * cos(2 * π * s0)

hh1 = Symbolics.derivative(hh, s0)
hh2 = Symbolics.derivative(hh, s1)
hh12 = Symbolics.derivative(hh1, s1)

function mybuild(x)
    aa = build_function(x, k, [s0, s1])
    eval(aa)
end

const h, h1, h2, h12 = mybuild.([hh, hh1, hh2, hh12])

function e(x0, x1, x2)
    h2(k, [x0, x1]) + h1(k, [x1, x2])
end

function getp(x...)
    u, up, du, dup = x
    1 / (h12(k, [u, up]) * du * dup)
end


function gen_solvem()
    a = [1 / (1 - exp(-2 * pi * ω * im * i)) for i in -n/2:-1]
    b = [1 / (1 - exp(-2 * pi * ω * im * i)) for i in 1:n/2-1]
    fftshift([a; [0]; b])
end


function gen_solvep()
    a = [1 / (exp(2 * pi * ω * im * i) - 1) for i in -n/2:-1]
    b = [1 / (exp(2 * pi * ω * im * i) - 1) for i in 1:n/2-1]
    fftshift([a; [0]; b])
end

gen_solvep()

##

const transcoeff1 = fftshift([exp(2 * pi * ω * im * i) for i in -n/2:(n/2)-1])
const transcoeff2 = fftshift([exp(-2 * pi * ω * im * i) for i in -n/2:(n/2)-1])
const dcoeff = fftshift([2 * pi * im * i for i in -n/2:(n/2)-1])
const dcoeffp = fftshift([2 * pi * im * i * exp(2 * pi * ω * im * i) for i in -n/2:(n/2)-1])
const solvem = gen_solvem()
const solvep = gen_solvep()
const id = [i / n + 0 * im for i in 0:n-1]
const P = plan_fft(id)
const iP = plan_ifft(id)


function pm(x, y, μ)
    c = similar(x)
    for i in eachindex(c)
        c[i] = x[i] * (y[i] + μ)
    end
    c
end

function iter!(u0v, u0f)
    u0pwf = transcoeff1 .* u0f
    u0pwf[1] = u0pwf[1] + n * ω
    u0pwv = iP * u0pwf
    u0mwf = transcoeff2 .* u0f
    u0mwf[1] = u0mwf[1] - n * ω
    u0mwv = iP * u0mwf
    du0f = dcoeff .* u0f
    du0v = iP * du0f
    du0pwf = dcoeffp .* u0f
    du0pwv = iP * du0pwf
    gv = -(1 .+ du0v) .* map(e, u0mwv + id, u0v + id, u0pwv + id)
    gf = P * gv
    ψv = iP * (gf .* solvem)
    pv = map(getp, id + u0v, id + u0pwv, 1 .+ du0v, 1 .+ du0pwv)
    μ = -sum(pv .* ψv) / sum(pv)
    wf = solvep .* (P * (pm(pv, ψv, μ)))
    wv = iP * (wf)
    preu0v = u0v + ((1 .+ du0v) .* wv)
    preu0f = P * (preu0v)
    for i in eachindex(preu0v)
        u0v[i] = preu0v[i]
        u0f[i] = preu0f[i]
    end
end

function runcode()
    u0v = [ω + 0im for i in 1:n]
    u0f = P * u0v
    for i in 1:100
        iter!(u0v, u0f)
    end
    u0v, u0f
end

@time aa = runcode();

直接法 Mathematica 代码:

n = 100;
\[Omega] = N[(Sqrt[5] - 1)/2];
k = 0.6;
h[{x0_, x1_}] = (1/2.) (x1 - x0)^2 - (k/(2. \[Pi])^2) Cos[
     2 \[Pi] x0];
e[um_, u_, up_] = D[h[{um, u}] + h[{u, up}], u];
h11[{x0_, x1_}] = D[h[{x0, x1}], x0, x0];
h22[{x0_, x1_}] = D[h[{x0, x1}], x1, x1];
h12[{x0_, x1_}] = D[h[{x0, x1}], x0, x1];
V = Table[
    ToExpression["v" <> ToString[i]], {i, 0, 
     n - 1}];(*未知函数V在区间[0,1]上的等分布*);
myfft[l_] := 
  Table[Sum[l[[i]]*E^(-2. \[Pi]*I*(i - 1)*(j - 1)/n), {i, 1, n}]/
    n, {j, 1, n}];
myifft[l_] := 
  Table[Sum[l[[i]]*E^(2. \[Pi]*I*(i - 1)*(j - 1)/n), {i, 1, n}], {j, 
    1, n}];
myfftshift[l_] := 
  Module[{a, b}, b = l[[n/2 + 1 ;; n]]; a = l[[1 ;; n/2]];
   Join[b, a]];
transcoeff1 = 
  myfftshift@
   Table[E^(2. \[Pi] I i \[Omega]), {i, -n/2, 
     n/2 - 1}];(*产生平移+\[Omega]的Fourier系数乘法列表,后N/2-1Subscript[个数遵循y, \
N]=Conjuate[Subscript[y, N-n]]的规则*)
transcoeff2 = 
  myfftshift@Table[E^(-2. \[Pi] I i \[Omega]), {i, -n/2, n/2 - 1}];
Vf = myfft[V];
Vpf = Times[Vf, transcoeff1];
Vpv = myifft[Vpf];(*利用离散Fourier变换求得v(\[Theta]+\[Omega])在[0,1]上的均匀分布值*)
Vmf = Times[Vf, transcoeff2];
Vmv = myifft[Vmf];
finaleq[umv_, uv_, upv_, Vm_, V_, 
   Vp_] := (h11[{uv, upv}] + h22[{umv, uv}]) V + h12[{uv, upv}]*Vp + 
   h12[{umv, uv}]*Vm;
SetOptions[Fourier, FourierParameters -> {-1, -1}];
SetOptions[InverseFourier, FourierParameters -> {-1, -1}];
geneq[eq_, v_] := Table[eq[[i]] == v[[i]], {i, 1, Length[eq]}];
id = Table[i/n, {i, 0., n - 1}];
u0v = \[Omega] & /@ id;

u0f = Fourier[u0v];(*利用离散Fourier变换求得v(\[Theta]-\[Omega])在[0,1]上的均匀分布值*)

Do[u0pf = 
  MapAt[# + \[Omega] &, MapThread[Times, {transcoeff1, u0f}], 
   1];(*考虑的映射写成\[Theta]+u(\[Theta])的形式,因此该映射的平移为\[Theta]+\[Omega]+u(\
\[Theta]+\[Omega]),其Fourier系数的第一个,即平均值需要加\[Omega]*)
 u0pv = InverseFourier[Take[u0pf]];
 u0mf = MapAt[# - \[Omega] &, MapThread[Times, {transcoeff2, u0f}], 1];
 (*考虑的映射写成\[Theta]+u(\[Theta])的形式,因此该映射的平移为\[Theta]-\[Omega]+u(\
\[Theta]-\[Omega]),其Fourier系数的第一个,即平均值需要减\[Omega]*)
 u0mv = InverseFourier[Take[u0mf]];
 de = MapThread[
   finaleq, {u0mv + id, u0v + id, u0pv + id, Vmv, V, Vpv}];
 
 ee = -MapThread[e, {u0mv + id, u0v + id, u0pv + id}];
 eq = geneq[de, ee];
 vv = V /. NSolve[eq, V][[1]];
 u0v = u0v + vv;
 u0f = Fourier[u0v];, 5](*迭代太多次会发散,在精确性以及稳定性上,远远不如使用同调方程迭代的方法*)

 twistmap[{x0_, y0_}] := {Mod[x0 + y0 + k/(2 \[Pi]) Sin[2 \[Pi] x0], 
    1], y0 + k/(2 \[Pi]) Sin[2 \[Pi] x0]};
getdatax0 = id + Re@u0v;
getdatay0pref = 
  MapAt[# + \[Omega] &, MapThread[Times, {transcoeff1, u0f}], 1];
getdatay0prev = 
  Re@InverseFourier[Take[getdatay0pref], 
    FourierParameters -> {-1, -1}];
getdatay0 = 
  MapThread[#2 - #1 - k/(2 \[Pi]) Sin[2 \[Pi] #1] &, {getdatax0, 
    id + getdatay0prev}];
data = Transpose[{Mod[getdatax0, 1], getdatay0}];
f = Interpolation[data];
{min, max} = f[[1, 1]];
Show[Plot[f[x], {x, min, max}, Frame -> True, AspectRatio -> 0.8, 
  PlotStyle -> {Thickness[0.0028], Black}, FrameLabel -> {x, y}, 
  FrameStyle -> 
   Directive[Black, FontSize -> 15, FontFamily -> "Times New Roman"]],
  ListPlot[NestList[twistmap, data[[1]], 100], PlotStyle -> {Red}]]