• Matlab:Toeplitz矩阵-向量乘法的快速傅里叶(FFT)算法


    一、$ t Toeplitz$矩阵与循环($ t Circulant$)矩阵

    定义

    为$n imes n$阶循环矩阵。

    定义 $T_n(i,j)=t_{j-i} $  为$n imes n$ 阶$ t Toeplitz$矩阵

    通过令矩阵$B_n=$

    从而构造出$2n imes 2n$阶循环矩阵

    假设有一$n imes 1$阶列向量$f u$

    其中,$C_{2n}$可以由快速傅里叶对角化

    其中$f c$表示$C_{2n}$矩阵的第一列元素,$f F$ 表示快速傅里叶($ t fft$)变换,$f F^{-1}$ 表示快速傅里叶($ t ifft$)逆变换。进一步可写成

    因此,计算$f T_n u$等价于计算

    查阅文献我们知道,直接计算$f T_n u$的存储量和计算量分别为$O(n^2)$和$O(n^3)$,但是利用快速傅里叶计算可以将存储量和计算量分别降为$O(n)$和$O(n log_2 n)$.

    从以下数据可以更直观的看出FFT显著的优点

     二、数值应用

    • 考虑一维椭圆方程

    $$-Delta u=f(x),qquad a<x<b, ag{1}$$

    满足齐次$Dirichlet$边界条件。

    对$xin [a,b]$一致网格剖分:$a=x_0<x_1,cdots,x_M=b$,$h=frac{b-a}{M}$。$U,u$分别表示数值解和真解。应用二阶中心差分逼近二阶导数

    $$Delta u(x_i)= frac{u(x_{i-1})-2u(x_i)+u(x_{i+1}) }{h^2}+O(h^2). ag{2}$$

    由(2)式可得求解方程(1)的数值格式的矩阵形式

    $$A{f U}=widehat{f}. ag{3}$$

    其中

    $$A= t -frac{1}{h^2}toeplitz([-2,1,zeros(1,M-3)]),$$

    $$widehat{f}=(   f(x_1),f(x_2),cdots,f(x_{M-1})   )^T.$$

    $${f U}=(u_1,u_2,cdots,u_{M-1})^T.$$

    • 考虑二维椭圆方程

    $$-Delta u=f({f x,y}),qquad {(f x,y)}in (a,b) imes (c,d), ag{4}$$

    对$xin [a,b]$一致网格剖分:$a=x_0<x_1,cdots,<x_{M_1}=b$,$h_1=frac{b-a}{M_1}$,$c=y_0<y_1,cdots,<y_{M_2}=d$,$h_2=frac{d-c}{M_2}$。$U,u$分别表示数值解和真解。应用二阶中心差分逼近二阶导数

    $$Delta u(x_i,y_j)= frac{u(x_{i-1},y_j)-2u(x_i,y_j)+u(x_{i+1},y_j) }{h_1^2}+ frac{u(x_i,y_{j-1})-2u(x_i,y_j)+u(x_i,y_{j+1}) }{h_2^2}+O(h_1^2+h_2^2). ag{5}$$

     由(5)式可得求解方程(4)的数值格式的矩阵形式

    $$A{f U}=widehat{f}. ag{6}$$

    其中

    $$A_1= t toeplitz([-2,1,zeros(1,M_1-3)]),$$

    $$A_2= t toeplitz([-2,1,zeros(1,M_2-3)]),$$

    $$ A_x = - t frac{1}{h_1^2} I_{M_2-1}  igotimes  A_1 ,mbox{(非toeplitz矩阵)}$$

    注意到:

    $$  I_{M_2-1}  igotimes  A_1U = reshapeBig( A_1 reshape( U,M_1-1,M_2-1 ),( M_1-1 )(M_2-1),1 Big). $$

    $$ A_y = - t frac{1}{h_2^2}  A_2 igotimes I_{M_1-1} ,$$

    $$A = A_x+A_y,$$

    $$widehat{f}=(   f(x_1,y_1),f(x_2,y_1),cdots,f(x_{M_1-1},y_1) , f(x_1,y_2),f(x_2,y_2),cdots,f(x_{M_1-1},y_2), cdotscdots, f(x_1,y_{M_2-1}),f(x_2,,y_{M_2-1}),cdots,f(x_{M_1-1},,y_{M_2-1}) )^T.$$

    $${f U}=(u_{1,1},u_{2,1},cdots,u_{M_1-1,1},u_{1,2},u_{2,2},cdots,u_{M_1-1,2},cdotscdots,u_{1,M_2-1},u_{2,M_2-1},cdots,u_{M_1-1,M_2-1})^T.$$

     由数值格式(3),(6)显然可知,当空间网格剖分很大时,矩阵乘法的计算量将会十分昂贵,因此利用FFT算法是很有必要的。接下来介绍一种有效的线性迭代法-共轭梯度法(CGS)

    三、数值例子

    • case $I$(1D) : 真解:

    $$ u = sin(x),qquad xin( 0,pi ), $$

    分别应用直接法和FFT方法的实验结果见下图

    • case $II$(2D) : 真解:

    $$ u = sin(x)sin(y),qquad (x,y)in( 0,pi )^2, $$

    分别应用直接法和FFT方法的实验结果见下图

     

    从数值实验结果可以直观的看出,FFT的计算效率是惊人的!

  • 相关阅读:
    STM32F0库函数初始化系列:进入STOP模式,外部中断唤醒
    STM32F0库函数初始化系列:ADC
    STM32F0库函数初始化系列:PWM输出
    VBS病毒实验
    AWVS漏洞扫描教程之扫描方式
    利用AWVS扫描Web漏洞
    命令执行漏洞靶场练习二
    命令执行漏洞靶场练习一
    CSRF POST型
    RainbowCrack彩虹表破解密码hash
  • 原文地址:https://www.cnblogs.com/xtu-hudongdong/p/11305766.html
Copyright © 2020-2023  润新知