LS文法与Julia逃逸时间算法改进

优化算法改进策略总结
优化算法改进策略总结的关键是根据具体问题的特点，选择合适的改进策略和技巧。

下面总结几种常见的优化算法改进策略：
1.贪心策略：贪心算法选择局部最优解，并希望通过不断选择
局部最优解来达到全局最优解。

贪心策略适用于那些具有贪心选择性质的问题。

2.动态规划：动态规划通过将原问题划分为多个子问题，并保
存子问题的解，通过递推求解子问题来得到原问题的解。

动态规划适用于具有重叠子问题和最优子结构的问题。

3.分支界定：分支界定通过建立一个解空间树，将搜索过程转
化为对解空间树的遍历，通过剪枝操作来减少搜索空间。

分支界定适用于具有可行解空间结构的问题。

4.回溯法：回溯法通过试探和回溯的方式来寻找问题的解，它
适用于具有多个可能解，并且每个可能解满足一定的约束条件的问题。

5.深度优先搜索：深度优先搜索通过不断地向前搜索到不能再
继续搜索为止，然后回退到上一个节点，再继续搜索。

深度优先搜索适用于解空间较大，但解的深度较小的问题。

6.广度优先搜索：广度优先搜索通过不断地将当前节点的所有
相邻节点入队，然后按照队列中的顺序进行遍历，直到找到目标节点或者遍历完所有节点。

广度优先搜索适用于解空间较小，
但解的广度较大的问题。

总的来说，对于优化算法的改进策略，需要根据具体问题的特点进行选择，针对问题的特点使用合适的算法和技巧，以提高算法的效率和准确性。

齐齐哈尔大学综合实践题目Julia分形及其Java编程学院理学院专业班级信科121班学生姓名指导教师成绩Julia 分形及其Java 编程由于本学期分形学老师所讲分形主要是以MATLAB 为例，所画出的分形图像效果并不理想，尤其经过几次放大后图像就会失真。

而Java 在分形的应用上效果比MATLAB 好很多，因此下文主要介绍Julia 分形图，实例是以Java 为基础。

分形是近几十年发展起来的一门新的数学分支，它涉及的领域非常之广，有物理学、数学、化学、生物学、医学、地震学、地貌学、冶金学、材料学、哲学、经济学、社会学等等.分形的出现正在改变科学家观察自然界的传统方式,目前已对当今数学乃至整个科学界产生了巨大的影响.本文主要对分形几何中的四元数进行研究，运用四元数绘制二维和三维Mandelbrot 集和Julia 集，并用Java 语言编程实现.先后介绍了分形的产生、Mandelbrot 集和Julia 集、四元数分形和用四元数绘制三维Mandelbrot 集和Julia 集的数学理论，关键词：分形，四元数，Julia 集，Java 程序设计一、Mandelbrot 集与Julia 集1、 Mandelbrot 集1980年, Mandelbrot 给世人提供了一幅无与伦比的杰作: Mandelbrot 集.其创作过程如下:令c z z f +=2)(, 其中C c z ∈,,z 是复变量，c 是复常数.对变换f 施行逃逸时间法, 得到如下迭代公式[1]:pq n n c z z +=+21 (1.1)式中)0,0(0=z ,pq c 为计算机荧屏位于),(q p 位置的象素. 于是(1.1)式成了pq pq pq pq pq n c c c c c z z ++++++=+22222201)))))((((( (1.2)给定N 为一个正整数, 比如等于255. 当象素位于),(q p 且N n =时,n z 仍然小于预设的一个阈值K , 则在),(q p 位置着色为1(蓝色), 否则当N n <时, 已有K z n ≥, 则在),(q p 位置描色为n . 如此),(q p 遍历整个荧屏后, 便画出了一幅Mandelbrot 集. 该集合的坐标如图2-1所示.图2-1 Mandelbrot 集的坐标图为什么说Mandelbrot 集是分形呢? 实在是它的层层嵌套中有很多很多的自相似部分.部分经逐级放大后, 又出现了一个Mandelbrot 湖.经典的M-集是由映射c z +2得来, 人们自然会采用更多的函数, 从而得到各种各样的M-集, 又称广义M-集.c z z n n +=+cos 1;c e z z c n n ++=+)arcsin(21; c c z z n n ++=+)cosh arctan(81.2、 Julia 集Gaston Julia(1893－1978)，法国数学家.1919年，他在第一次世界大战时受了伤，住院期间他潜心研究了迭代保角变换c z z n n +=+21. 这种复平面上的变换能生出一系列令人眼花缭乱的图形变化. 当时没有电子计算机，不能像现在那样把如此美妙绝论的图案奉献于世. 因此他的工作并不为世人重视.虽然产生Julia 集和Mandelbrot 集的变换都是3,2,1,0,21=+=+n c z z n n (1.3)但这里的常数c 却是任意复数, 变元0z 是计算机荧屏上的每个象素. 当0z 遍历象素),(q p 的所有点且对公式(1.3)运用逃逸时间法后, 便得到一幅Julia 集c J 了.令c= -0.65175; 0.41850;便得到 图2.2图2-2Julia 集(c= -0.65175; 0.41850;)由于Mandelbrot 集(有时简称其为M-集)和Julia 集都源于同一个变换, 因此它们之间必定有非常复杂的关系. 由于每一个常数 c 都对应一个c J , 而M-集上的每个点都是一个c , 所以M-集合的所有点就对应着数以万计的c J . 图2-3正显示了它们两者之间的这种关系. 看得出, 相近的c 值, 对应的c J 也就较为相似.图2-3 Mandelbrot 集和Julia 集从图2.3可以看出, 凡是M-集的边界点, 其对应的c J 就显有分枝状. 这里面有太多的奥妙.由于一幅Julia 集完全依赖于常数c , 所以我们常常把它简记为c J .图2-4 c 不同的Julia 集与M-集一样, Julia 集也有其广义集. 图2.4的四幅图分别对应于函数:1. )sin(cos sin 2z z c z =2. )cos(πz c z ⋅=3. πc z z z +⋅=cos )log(cosh4. πc z z z +⋅=cos )tanh(cos.二维Julia 集，用牛顿迭代算法方法进行图形绘制，绘制范围：5.1～5.1:-x 0.2～0.2:-yC= (0.67,0.48)绘制结果如图2-5所示：图2-5 二维Julia 集之二二、四元数分形与四元数Mandelbrot 集和Julia 集1、 四元数基本理论我们都熟悉平面上复数z , 其中bi a z +=, 而1-=i . 无疑, 十八世纪以前创立的复数是数学史上的一件大事. 那末是否有高维的复数呢? 所谓的超复数(hypercomplex 或supercomplex)是这样定义的:令ck bj ai w z +++=, 其中i 和k j ,都是虚数, 它们满足下述运算要求:11=-====-=-==-=-==ijk kk jj ii j ik i kj k ji j ki i jk k ij看得出, 它们的乘法满足交换律.两个超复数的乘法公式是:令 k z j y i x w h 11111+++=和k z j y i x w h 22222+++=.则+---=)(2121212121z z y y x x w w h h +--+i z y y z x w w x )(21212121 +-+-j z x y w x z w y )(21212121 kz w y x x y w z )(21212121+++其它运算法则就不再赘述.1843年, 爱尔兰数学家William R. Hamilton(哈密尔顿)发明了四元数(quaternion). 一个四元数q 的定义是这样的:zk yj xi w q +++= (2.1)其中z y x w ,,,是实数,k j i ,,是虚数, 且有1222-====ijk k j i (2.2)q 的模2222z y x w q +++=. (2.3)一个四元数Q 可以由平面上的两个复数v u ,来表示:⎥⎦⎤⎢⎣⎡-+-++=⎥⎦⎤⎢⎣⎡-=bi a di c di c bi a u v v u Q (2.4) 其中d c b a ,,,是实数.v u ,分别是复数v u ,的共轭.一个四元数也能用[]w z y x q ,,,= (2.5)来表示.一个平面上的复数由实部和虚部构成: bi a z +⋅=1, 一个四元数Q 同样也能由若干部分线性组合而成:zK yJ xI wU Q +++= (2.6)其中⎥⎦⎤⎢⎣⎡≡1001U (2.7)⎥⎦⎤⎢⎣⎡-≡i i I 00 (2.8) ⎥⎦⎤⎢⎣⎡-≡0110J (2.9) ⎥⎦⎤⎢⎣⎡≡00i i K (2.10) 于是U K J I -===222 (2.11)也就是说,K J I ,,是矩阵方程 U X -=2的解, 是负单位矩阵的平方根.一个四元数整系数基的线性组合也叫Hamilton 整数. 在4R 空间, 四元数的基是如下四个:⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=10000010001000011 (2.12) ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=0100100000010010i (2.13) ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=0001001001001000j (2.14) ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=0010000110000100k (2.15) 与超复数不同, 四元数的三个虚数之间的运算并不遵从乘法交换律, 其运算规则如下:k ji ij =-= (2.16) i kj j =-= (2.17) j ik ki =-= (2.18)看起来很象三维空间直角坐标系中三个单位向量k j i ,,的叉乘关系.设 zk yj xi w q +++=, 则其四元共轭为zk yj xi w q ---=.其加法遵从一般规律:k z z j y y i x x w w q q )()()()(2121212121+++++++=+.设],,,[],,,,[22221111z y x w q z y x w p ==, 其乘法服从)(21212121z z y y x x w w pq ---= i y z z y w x x w )(21212121-+++ j x z w y z x y w )(21212121++-+ kw z x y y x z w )(21212121+-++q 的模仍然遵从一般复数关系.)(q q q q q norm == (2.19)且等于公式(2.3).一个四元数可以写成一个数量加上一个向量).,(v w q=其中向量 [,,]v x y z =. 如此一来, 两个四元数的乘法就变得较为简单:),(),(221121v w v w q q⋅=),(2112212121v v v w v w v v w w⨯++⋅-= 四元数的除法也遵从复数关系qq qp q p p p p p ==-和1 从几何上来讲, 四元数代表着时间加三维空间. 如果固定实数w 为常数, 则这个四元数就是三维空间的一个变量.0)(lim=→qq f q2 四元数Mandelbrot 集和Julia 集记q 和c 为四元数, 按第二章中Mandelbrot 集的迭代法,2,1,0,21=+=+n c q q n n (2.20)就得到一幅3维四元数(把实数w 看作常数, 并认定其为第四维)Mandelbrot 集了.四元数Julia 集是4维空间中的Julia 集，同样使用公式,2,1,0,21=+=+n c q q n n其中),2,1,0( =n q n 和c 都是四元数. 生成四元数Julia 集的方法与普通复数情况类似，只需迭代公式(2.20), 只是c 是一个固定的四元数而已. 迭代时要观察 ||n q 的敛散情况.把第四维看作时间，可以认为我们生成的是一个三维动画片. 实际操作时,一般取第四维为常数，这样我们得到的是静止的三维图像. 程序主体结构与二维情况也类似，我们用多重循环扫描一个三维立方体内部所有的点，就知道该点是否属于Julia 集了. 做法如下: 对每一个y对每一个 x 对每一个z作迭代,2,1,0,21=+=+n c q q n n其中],,,0[n n n n z y x q =. 当M n =时, n q 仍属于集合的话就退出这一层循环我们并不计算z 轴方向上的所有点，那样速度太慢. 对于给定的一组),(y x ，只要沿z 轴方向找到第一个属于Julia 集合的点就够了, 其它的点被挡住，让看不见. 这样，就大大提高程序的运行速度. 也有广义的四元数Julia 集. 不过由于以四元数为自变量的函数计算过于困难, 目前还仅仅处理q 的多项式.结 论在分形几何中，许多重要的分形是由迭代产生的.因为迭代可以使一些看似简单的函数产生惊人的复杂性，Julia 集就是其中一种.由于迭代函数的多样性，Julia 集可以在计算机图形工具的辅助下呈现为色彩斑斓、结构优美的分形图案，因此可广泛应用于纺织印染、广告设计、服装设计、装潢设计以及计算机美术教学等领域.可见，研究Julia 集的生成算法具有重要的理论意义和实际价值.逃逸时间算法是生成Julia 集的经典算法，它具有基本原理简单、绘图精度高、占用内存少等优点附：代码。

ls自适应算法原理自适应算法（Least Squares Adaptation, LS算法）是一种用于实时信号处理和自适应滤波的方法。

它通过对输入信号和期望输出信号之间的误差进行最小二乘拟合，来逐步调整滤波器的系数，从而达到滤波的目的。

本文将详细介绍ls自适应算法的原理及其应用。

1. 概述自适应滤波器是一种根据输入信号的统计特性自动调整其传递函数的滤波器。

相比于传统的滤波器，自适应滤波器可以更好地适应不断变化的信号环境，提供更好的性能。

LS算法是一种常用的自适应滤波算法，其核心思想是通过最小化误差平方和来不断调整滤波器的系数。

2. LS算法原理LS算法通过最小二乘法来估计滤波器的系数。

假设有一个长度为N的滤波器系数向量W，输入信号向量X以及期望输出信号向量D。

LS算法的目标是使得实际输出信号Y与期望输出信号D之间的误差最小，即最小化损失函数L：L = ||D - Y||²其中，||.||表示向量的范数。

为了最小化损失函数，我们要求导损失函数L关于滤波器系数向量W的偏导为零：∂L / ∂W = 0对上式求导得到：∂L / ∂W = -2X^\mathrm{T}(D - XW)其中，^\mathrm{T}表示矩阵的转置。

令上式等于零，可以得到滤波器系数的估计值W：Ŵ = (X^\mathrm{T}X)⁻¹X^\mathrm{T}D3. LS算法应用LS算法在自适应滤波中有广泛的应用，例如回声消除、自适应降噪和信号预测等领域。

下面以自适应降噪为例，介绍LS算法的应用过程。

（1）收集数据首先，需要收集包含有噪声的输入信号向量X和对应的期望输出信号向量D。

这可以通过传感器获取信号数据，并根据信号的统计特性进行采样和处理得到。

（2）初始化滤波器系数为了使用LS算法进行自适应降噪，需要初始化滤波器系数向量W。

一种常用的初始化方法是将滤波器系数置为零向量。

（3）使用LS算法调整滤波器系数根据LS算法的原理，在收集到的数据上应用LS算法，通过计算得到滤波器系数估计值Ŵ。

FFT算法的的改进FFT(Fast Fourier Transform)是一种高效的算法，用于计算离散傅里叶变换(DFT)。

它通过将DFT的时间复杂度从O(N^2)降低到O(NlogN)，在信号处理、图像处理和数值计算等领域得到广泛应用。

然而，有许多对FFT的改进方法，旨在进一步提高其性能和功能。

在本文中，将介绍一些常见的FFT算法改进。

1. Cooley-Tukey算法Cooley-Tukey算法是最早提出的FFT算法之一，它将DFT的计算过程分解成更小规模的DFT计算。

该算法利用了傅里叶变换的周期性质，并将复杂度从O(N^2)降低到O(NlogN)。

Cooley-Tukey算法的关键思想是将DFT分解为两个规模较小的DFT。

具体而言，该算法将输入序列分为奇数索引和偶数索引，然后分别对它们进行DFT计算。

接下来，再将结果通过对称性质合并起来。

该算法的递归实现可以进一步降低计算复杂度。

2.快速数论变换(NTT)快速数论变换(NTT)是对Cooley-Tukey算法的拓展，它应用于模数为p的多项式乘法问题。

NTT利用了复杂指数的周期性质，并适用于计算形式化指数有规律的DFT。

与FFT类似，NTT的核心思想是将DFT分解为较小规模的DFT计算。

不过，NTT使用了不同的运算规则，因此需要适应于特定的模数p。

通过结合合适的数论性质，NTT可以在O(NlogN)的时间复杂度内计算出DFT。

3. Radix-2算法Radix-2算法是基于Cooley-Tukey算法的改进版本，它在进行DFT计算时要求输入序列的长度必须是2的幂次。

Radix-2算法利用了二进制的特性，将DFT计算分解为若干层级的蝶形运算。

每一层级都将输入序列分为两部分，并在频域上进行乘法和加法运算。

通过逐层级地进行运算，Radix-2算法可以快速计算出DFT。

4. 快速Hadamard变换(FHT)快速Hadamard变换(FHT)是对FFT的另一种改进算法，它适用于实数域上的信号处理。

经典Julia 集和Mandelbrot 集的几何画板画法重庆市万州二中 向忠复平面内由二次函数2()f z z c =+构造的迭代点列 {}n z ：10(),,1,2,3,n n z f z z z n −===能够生成分形图形。

当复数c 为定值(|c|<2，c ≠0)时，集合0{|{}}c n J z z =有界称为Julia 集；而对于复数z 为定值(通常z=0)时，集合{|{}}z n M c z =有界称为Mandelbrot 集。

在绘制Julia 集或Mandelbrot 集时，对于集合c J 或z M 外的元素，复数()n f z 的模可能趋于∞，这样在复平面相应位置处，图形将会出现空洞。

为使图形铺满整个画面，使分形图形更加清晰漂亮，我们常常采用逃逸时间算法作分形图。

设置一个区域A(通常设为圆|z|≤2)，对于充分大的整数N ，当点0z 经过n+1次(n ≤N)迭代后的1n z +跳出了区域A(即1||2n z +>，我们就说点1n z +逃逸了)，记录et 为n 的值；而如果点z 经过N 次迭代后的点仍未超出这个区域，就记录et 为N ，这个值et 就叫逃逸时间，而em=||et z 叫逃逸半径。

利用逃逸时间和逃逸半径关联着色参数对点0z 或c 着色，即能描绘出集合c J 或z M 的图形。

这就是逃逸时间算法的基本思想。

通俗的讲，经典Mandelbrot 集和Julia 集，就是将复平面内的圆|z|=2，通过逃逸时间算法下的迭代2z z +c →生成的复杂图形。

几何画板制作经典Mandelbrot 集和Julia 集的步骤：1.设置角度为弧度,作点z 、c 的横纵坐标：新建参数xz 、yz 、xc 、yc ；计算z 的模方||z||= xz*xz+yz*yz ；计算阈判断真值p=sgn(1-sgn(||z||-4))；2.计算点z 在迭代2z z +c →下的象点Z 的横纵坐标：xZ=xz*xz-yz*yz+xc 、yZ=2xz*yz+yc；3.构建点z的逃逸变换点Z'，即若点z在逃逸阈内，则p=1，点z变换为点Z'；否则p=0，点z停止于z：计算Z'的横纵坐标：xZ'=xz+p*(xZ-xz)、yZ'=yz+p*(yZ-yz)；4.为记录逃逸时间和逃逸半径，新建参数t0=0，计算t0+p；作点E(t0, ||z||)；5.新建n=3，作{xz、yz、t0}→{xZ'、yZ'、t0+p}的深度为n的迭代，继而作点E的迭代象终点iE，并度量其横纵坐标et(逃逸时间)和em(逃逸半径)；6.计算RGB着色参数：计算s=.05(et-log(.5abs(ln(em))))，R=sin(s)，G=sin(3s)，B=cos(2s)；7.作点pixel，度量其横纵坐标x pixel、y pixel；对点pixel进行RGB着色并作颜色变换。

高阶Julia集的逃逸时间算法的计算机实现作者：董洁冯毅宏刘冬莉来源：《硅谷》2011年第18期摘要： Julia集是分形理论中具有重要地位的集合，近些年来有很多对于变换f（z）=z2+c 生成分形图形的扩展研究，主要针对高阶非线性复映射迭代函数f（z）=zn+c，给出利用具体的逃逸时间算法生成分形图的具体步骤以及生成的Julia集图案。

关键词： Julia集；逃逸时间算法；分形；函数迭代中图分类号：TP391.41 文献标识码：A 文章编号：1671－7597（2011）0920183－011 概述分形几何中，许多重要的分形是由迭代产生的。

因为迭代可以使一些看似简单的函数产生色彩斑斓、千变万化而又具有任意高分辨率结构的艺术图案，Julia集即为其中一种。

近些年来有很多对于变换f（z）=z2+c生成分形图形的扩展研究，取得了很多理论成果和美丽的图形，但提到实现的算法时却步骤简单，本文即是在对迭代公式的扩展：f（z）=zn+c（n>2）的基础上，给出了实现算法的详细步骤。

2 Julia集概念与性质复平面C上的简单变换（函数）可以生成一些结构复杂的集合。

Julia集提供了利用简单函数的迭代过程却能生成复杂结构的集合的典型例子。

其定义如下：设是阶数大于1的多项式，表C中那些轨道不趋于无穷点的点的集合，即：，称此集为对应于f的充满的Julia集，的边界称为多项式f的Julia集，记为即：。

[1]在计算机上绘制Julia集图形，需将复平面上的数据引入到平面直角坐标系中。

对于高次复迭代函数f（z）=zn+c（n＞2），迭代函数形式可表示为：zk+1=zkn+c（zk，c∈C，k=0，1，2，…）。

其中，Zk为第k次迭代后的复数xk+yki，C为常数cx+cyi，则每次给定的初始值Z0，代入公式可得到Z1，经过k次迭代后可得到复数序列：Z2，Z3，……，Zk，当c≠0且其实部和虚部都取较小的实数时，这个序列将有如下情况：①如果Z0的模也很小，当迭代次数达到足够多的时候，则迭代序列接近吸引子。

基于Julia的科学计算算法优化与性能调试Julia是一种高性能动态编程语言，专为科学计算而设计。

它具有与C语言相媲美的性能，同时拥有易于使用的语法和丰富的库函数，使得它成为科学计算领域的瑰宝。

在进行科学计算时，编写高效的算法并对其性能进行调试优化是至关重要的。

本文将介绍如何基于Julia 进行科学计算算法的优化与性能调试。

1. 算法优化在进行科学计算时，编写高效的算法是提高计算效率的关键。

下面是一些优化算法的常用技巧：1.1 向量化操作在Julia中，向量化操作可以显著提高代码的执行速度。

通过使用向量和矩阵运算，可以避免使用循环，从而减少计算时间。

示例代码star：编程语言：julia# 未向量化操作result = zeros(n)for i in 1:nresult[i] = x[i] + y[i]end# 向量化操作result = x + y示例代码end1.2 避免不必要的内存分配在编写代码时，尽量避免不必要的内存分配操作。

可以通过预分配内存空间或者使用in-place操作来减少内存开销。

示例代码star：编程语言：julia# 避免不必要的内存分配result = zeros(n)for i in 1:nresult[i] = x[i] + y[i]end# 预分配内存空间result = similar(x)result .= x .+ y示例代码end1.3 使用原生Julia函数尽量使用原生Julia函数而不是调用外部库函数，因为原生Julia函数通常能够更好地利用Julia语言的特性和优化。

示例代码star：编程语言：julia# 使用原生Julia函数function mysum(x)s = zero(eltype(x))for xi in xs += xiendreturn send示例代码end2. 性能调试除了优化算法外，对代码的性能进行调试也是非常重要的。

Julia提供了一些工具和技巧来帮助我们进行性能调试。

julia解方程组摘要：一、Julia编程语言简介二、解方程组的基本方法1.高斯消元法2.矩阵求逆法3.迭代法三、Julia解方程组示例四、Julia代码实战五、总结与展望正文：Julia是一种高性能、通用、动态类型的编程语言，以其简洁的语法和强大的数值计算能力而著称。

在本文中，我们将介绍如何使用Julia解方程组。

首先，我们来了解一下Julia编程语言的基本概念。

Julia源于Python，但相较于Python，Julia在性能上有显著优势，特别是在科学计算方面。

Julia拥有丰富的库，可以轻松完成各种数学计算、图形绘制和数据分析等任务。

接下来，我们讨论解方程组的基本方法。

在Julia中，我们可以使用以下三种方法解方程组：1.高斯消元法：这是一种经典的数值解法，通过逐步消元将方程组转化为单变量方程，然后求解各个变量。

在Julia中，可以使用内置的` solve`函数实现高斯消元法。

2.矩阵求逆法：当方程组的系数矩阵满足一定条件时，可以采用矩阵求逆法求解。

在Julia中，可以使用`inv`函数求矩阵的逆。

3.迭代法：迭代法是一种简单且实用的解法，通过不断更新变量值直至收敛。

在Julia中，可以使用`while`循环实现迭代法。

接下来，我们通过一个示例来展示如何使用Julia解方程组。

假设我们有一个三元一次方程组：```3x + 2y + z = 12x - 4y + 2z = 10-x + 2y - z = -6```我们可以使用高斯消元法求解此方程组。

首先，设方程组为：```Ax = b```其中，A为系数矩阵，x为变量向量，b为常数向量。

我们可以通过以下步骤求解：1.构建系数矩阵A和常数向量b。

2.调用`solve`函数，求解变量向量x。

3.输出解x。

以下是Julia代码实战：```juliafunction solve_system_of_equations(A, b)n = size(A, 1)x = Array(Float64, n)for i in 1:nx[i] = b[i] / sum(A[i, i:n])endreturn xendA = [3, 1; 1, 2; -1, 2]b = [12; 10; -6]x = solve_system_of_equations(A, b)println(x)```运行上述代码，我们可以得到方程组的解：`[2.0, 1.0, 1.0]`。

分形几何中一些经典图形的Matlab画法（1）Koch曲线程序koch.mfunction koch(a1,b1,a2,b2,n)%koch(0,0,9,0,3)%a1,b1,a2,b2为初始线段两端点坐标,n为迭代次数a1=0;b1=0;a2=9;b2=0;n=3;%第i-1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中[A,B]=sub_koch1(a1,b1,a2,b2);for i=1:nfor j=1:length(A)/5;w=sub_koch2(A(1+5*(j-1):5*j),B(1+5*(j-1):5*j));for k=1:4[AA(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5),BB(5*4*(j-1)+5*(k-1)+1:5*4*(j-1)+5*(k-1)+5)] =sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4));endendA=AA;B=BB;endplot(A,B)hold onaxis equal%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别存%储在数组A,B中function [A,B]=sub_koch1(ax,ay,bx,by)cx=ax+(bx-ax)/3;cy=ay+(by-ay)/3;ex=bx-(bx-ax)/3;ey=by-(by-ay)/3;L=sqrt((ex-cx).^2+(ey-cy).^2);alpha=atan((ey-cy)./(ex-cx));if (ex-cx)<0alpha=alpha+pi;enddx=cx+cos(alpha+pi/3)*L;dy=cy+sin(alpha+pi/3)*L;A=[ax,cx,dx,ex,bx];B=[ay,cy,dy,ey,by];%把由函数sub_koch1生成的五点横、纵坐标A,B顺次划分为四组,分别对应四条折线段中%每条线段两端点的坐标,并依次分别存储在4*4阶矩阵k中,k中第i(i=1,2,3,4)行数字代表第%i条线段两端点的坐标function w=sub_koch2(A,B)a11=A(1);b11=B(1);a12=A(2);b12=B(2);a21=A(2);b21=B(2);a22=A(3);b22=B(3);a31=A(3);b31=B(3);a32=A(4);b32=B(4);a41=A(4);b41=B(4);a42=A(5);b42=B(5);w=[a11,b11,a12,b12;a21,b21,a22,b22;a31,b31,a32,b32;a41,b41,a42,b42];图1 V on Koch曲线（2）Levy 曲线程序levy.mfunction levy(n)% levy(16),n为levy曲线迭代次数%x1,y1,x2,y2为初始线段两端点坐标,nn为迭代次数n=16;x1=0;y1=0;x2=1;y2=0;%第i-1次迭代时由各条线段产生的新两条线段的三端点横、纵坐标存储在数组X、Y中[X,Y]=levy1(x1,y1,x2,y2);for i=1:nfor j=1:length(X)/3w=levy2(X(1+3*(j-1):3*j),Y(1+3*(j-1):3*j));[XX(3*2*(j-1)+1:3*2*(j-1)+3),YY(3*2*(j-1)+1:3*2*(j-1)+3)]=levy1(w(1,1),w(1,2),w(1,3),w(1,4) );[XX(3*2*(j-1)+3+1:3*2*(j-1)+3+3),YY(3*2*(j-1)+3+1:3*2*(j-1)+3+3)]=levy1(w(2,1),w(2,2),w( 2,3),w(2,4));endX=XX;Y=YY;endplot(X,Y)hold onaxis equal%由以(x1,y1),(x2,y2)为端点的线段生成新的中间点坐标并把(x1,y1),(x2,y2)连同新点横、纵坐%标依次分别存储在数组X,Y中function [X,Y]=levy1(x1,y1,x2,y2)x3=1/2*(x1+x2+y1-y2);y3=1/2*(-x1+x2+y1+y2);X=[x1,x3,x2];Y=[y1,y3,y2];%把由函数levy1生成的三点横、纵坐标X,Y顺次划分为两组,分别对应两条折线段中每条线%段两端点的坐标,并依次分别存储在2*4阶矩阵w中,w中第i(i=1,2)行数字代表第i条线段%两端点的坐标function w=levy2(X,Y)a11=X(1);b11=Y(1);a12=X(2);b12=Y(2);a21=X(2);b21=Y(2);a22=X(3);b22=Y(3);w=[a11,b11,a12,b12;a21,b21,a22,b22];图2 Levy 曲线（3）分形树程序tree.hfunction tree(n,a,b)% tree(8,pi/8,pi/8),n为分形树迭代次数%a,b为分枝与竖直方向夹角%x1,y1,x2,y2为初始线段两端点坐标,nn为迭代次数n=8;a=pi/8;b=pi/8;x1=0;y1=0;x2=0;y2=1;plot([x1,x2],[y1,y2])hold on[X,Y]=tree1(x1,y1,x2,y2,a,b);hold onW=tree2(X,Y);w1=W(:,1:4);w2=W(:,5:8);% w为2^k*4维矩阵,存储第k次迭代产生的分枝两端点的坐标, % w的第i(i=1,2,…,2^k)行数字对应第i个分枝两端点的坐标w=[w1;w2];for k=1:nfor i=1:2^k[X,Y]=tree1(w(i,1),w(i,2),w(i,3),w(i,4),a,b);W(i,:)=tree2(X,Y);endw1=W(:,1:4);w2=W(:,5:8);w=[w1;w2];end%由每个分枝两端点坐标(x1,y1),(x2,y2)产生两新点的坐标(x3,y3),(x4,y4),画两分枝图形,并把%(x2,y2)连同新点横、纵坐标分别存储在数组X,Y中function [X,Y]=tree1(x1,y1,x2,y2,a,b)L=sqrt((x2-x1)^2+(y2-y1)^2);if (x2-x1)==0a=pi/2;else if (x2-x1)<0a=pi+atan((y2-y1)/(x2-x1));elsea=atan((y2-y1)/(x2-x1));endendx3=x2+L*2/3*cos(a+b);y3=y2+L*2/3*sin(a+b);x4=x2+L*2/3*cos(a-b);y4=y2+L*2/3*sin(a-b);a=[x3,x2,x4];b=[y3,y2,y4];plot(a,b)axis equalhold onX=[x2,x3,x4];Y=[y2,y3,y4];%把由函数tree1生成的X,Y顺次划分为两组,分别对应两分枝两个端点的坐标,并存储在一维%数组w中function w=tree2(X,Y)a1=X(1);b1=Y(1);a2=X(2);b2=Y(2);a3=X(1);b3=Y(1);a4=X(3);b4=Y(3);w=[a1,b1,a2,b2,a3,b3,a4,b4];图3 分形树（4）IFS算法画Sierpinski三角形程序sierpinski_ifs.hfunction sierpinski_ifs(n,w1,w2,w3)%sierpinski_ifs(10000,1/3,1/3,1/3)%w1,w2,w3出现频率n=10000;w1=1/3;w2=1/3;w3=1/3;M1=[0.5 0 0 0 0.5 0];M2=[0.5 0 0.5 0 0.5 0];M3=[0.5 0 0.25 0 0.5 0.5];x=0;y=0;% r为[0,1]区间内产生的n维随机数组r=rand(1,n);B=zeros(2,n);k=1;% 当0<r(i)<1/3时,进行M1对应的压缩映射;% 当1/3=<r(i)<2/3时,进行M2对应的压缩映射;% 当2/3=<r(i)<1时,进行M3对应的压缩映射;for i=1:nif r(i)<w1a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);else if r(i)<w1+w2a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);else if r(i)<w1+w2+w3a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);endendendx=a*x+b*y+e;y=c*x+d*y+f;B(1,k)=x;B(2,k)=y;k=k+1;endplot(B(1,:),B(2,:),'.','markersize',0.1)图4 Sierpinski三角形（5）IFS算法画Julia集程序julia_ifs.hfunction julia_ifs(n,cx,cy)% julia_ifs(100000,-0.77,0.08)% f(z)=z^2+c,cx=real(c);cy=image(c);n=10000;cx=-0.77;cy=0.08;% z^2+c=z0,x=real(z0);y=image(z0);x=1;y=1;B=zeros(2,n);k=1;% A为产生的服从标准正态分布的n维随机数组A=randn(1,n);for i=1:nwx=x-cx;wy=y-cy;if wx>0alpha=atan(wy/wx);endif wx<0alpha=pi+atan(wy/wx);endif wx==0alpha=pi/2;endalpha=alpha/2;r=sqrt(wx^2+wy^2);if A(i)<0r=-sqrt(r);elser=sqrt(r);endx=r*cos(alpha);y=r*sin(alpha);B(1,k)=x;B(2,k)=y;k=k+1;endplot(B(1,:),B(2,:),'.','markersize',0.1)图5 Julia 集(6)逃逸时间算法画Sierpinski垫片程序sierpinski.hfunction sierpinski(a,b,c,d,n,m,r)%sierpinski(0,0,1,1,12,200,200)%(a,b),(c,d)收敛区域左上角和右下角坐标,m为分辨率% n为逃逸时间,需要反复试探,r逃逸半径a=0;b=0;c=1;d=1;n=12;m=200;r=200;B=zeros(2,m*m);w=1;for i=1:mx0=a+(c-a)*(i-1)/m;for j=1:my0=b+(d-b)*(j-1)/m;x=x0;y=y0;for k=1:nif y>0.5x=2*x;y=2*y-1;else if x>=0.5x=2*x-1;y=2*y;elsex=2*x;y=2*y;endif x^2+y^2>rbreak;endendif k==nB(1,w)=i;B(2,w)=j;w=w+1;endendendplot(B(1,:),B(2,:),'.','markersize',0.1)图6 Sierpinski三角形垫片(7)元胞自动机算法画Sierpinski三角形程序一维元胞自动机sierpinski_ca1.hfunction sierpinski_ca1(m,n)%sierpinski_ca1(1000,3000)m=1000;n=3000;x=1;y=1;t=1;w=zeros(2,m*n);s=zeros(m,n);s(1,fix(n/3))=1;for i=1:m-1for j=2:n-1if (s(i,j-1)==1&s(i,j)==0&s(i,j+1)==0)|(s(i,j-1)==0&s(i,j)==0&s(i,j+1)==1) s(i+1,j)=1;w(1,t)=x+3+3*j;w(2,t)=y+5*i;t=t+1;endendendplot(w(1,:),w(2,:),'.','markersize',1)图7.1 一维元胞自动机画Sierpinski三角形二维元胞自动机sierpinski_ca2.hfunction sierpinski_ca2(m,n)%sierpinski_ca2(400,400)m=400;n=400;t=1;w=zeros(2,m*n);s=zeros(m,n);s(m/2,n/2)=1;for i=[m/2:-1:2,m/2:m-1]for j=[n/2:-1:2,n/2:n-1]ifmod(s(i-1,j-1)+s(i,j-1)+s(i+1,j-1)+s(i-1,j)+s(i+1,j)+s(i-1,j+1)+s(i,j+1)+s(i+1,j+1),2)==1 s(i,j)=1;w(1,t)=i;w(2,t)=j;t=t+1;endendendplot(w(1,:),w(2,:),'.','markersize',0.1)图7.2 二维元胞自动机画Sierpinski三角形(8)IFS算法画Helix曲线程序helix_ifs.hfunction helix_ifs(n,w1,w2,w3)%helix_ifs(20000,0.9,0.05,0.05)%w1,w2,w3为出现频率n=20000;w1=0.9;w2=0.05;w3=0.05;M1=[0.787879 -0.424242 1.758647 0.242424 0.859848 1.408065];M2=[-0.121212 0.257576 -6.721654 0.05303 0.05303 1.377236];M3=[0.181818 -0.136364 6.086107 0.090909 0.181818 1.568035];x=0;y=0;% r为[0,1]区间内产生的n维随机数组r=rand(1,n);B=zeros(2,n);k=1;% 当0<r(i)<1/3时,进行M1对应的压缩映射;% 当1/3=<r(i)<2/3时,进行M2对应的压缩映射;% 当2/3=<r(i)<1时,进行M3对应的压缩映射;for i=1:nif r(i)<w1a=M1(1);b=M1(2);e=M1(3);c=M1(4);d=M1(5);f=M1(6);else if r(i)<w1+w2a=M2(1);b=M2(2);e=M2(3);c=M2(4);d=M2(5);f=M2(6);else if r(i)<w1+w2+w3a=M3(1);b=M3(2);e=M3(3);c=M3(4);d=M3(5);f=M3(6);endendendx=a*x+b*y+e;y=c*x+d*y+f;B(1,k)=x;B(2,k)=y;k=k+1;endplot(B(1,:),B(2,:),'.','markersize',0.1)图8 Helix曲线。

动态网络最短路径射线追踪算法中向后追踪方法的改进摘要：动态网络最短路径射线追踪算法中的向后追踪方法能够解决线性走时插值算法（LTI）向后追踪过程不稳定的问题，但是其计算效率较低.综合利用节点次级源的位置信息以及波的传播规律，提出了改进方法，排除了动态网络最短路径射线追踪算法向后追踪过程中存在的大量冗余计算.数值算例表明，改进的向后追踪方法具有较高的计算效率，是动态网络最短路径射线追踪算法中向后追踪方法的几倍至几十倍；若将改进后的向后追踪方法应用于动态网络最短路径射线追踪改进算法，则该算法的计算效率将提高一倍左右.关键词：射线追踪；线性走时插值；向后追踪方法；计算效率；初至波射线追踪中图分类号：P631 文献标识码：A文章编号：1674-2974（2016）05-0106-07Abstract：The backward tracing method of the shortest path ray tracing algorithm with dynamic networks can solve the unstability problem in the backward tracing procedure of theLTI （Linear Travel-time Interpolation）algorithm，but the computational efficiency of the method is low. This studypresented an improved method on backward tracing. According to the location information of the secondary sources for the nodes and the law of wave propagation，a large number of redundancy calculation are excluded in the backward tracing of the dynamic networks tracing algorithm. The numerical examples show that the improved method exhibits the higher computational efficiency. The calculation efficiency of the improved method is several times that of the backward tracing method of the dynamic networks tracing algorithm. When the improved method is applied to the improved algorithm of the shortest path ray tracing with dynamic networks，the computational efficiency of the algorithm can be increased by about 100 %.Key words：ray tracing；linear traveltime interpolation；improved algorithm；backward tracing；computational efficiency；first arrival ray tracing射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用.目前常用射线追踪方法主要有两点射线追踪算法（包括试射法以及弯曲法）[[1-3]、有限差分解程函方程法[[4-6]、最短路径法[[7-10]以及LTI（Linear Travel-time Interpolation）射线追踪算法[[11-22]等.其中，LTI 射线追踪算法因其计算精度较高、计算速度较快且适用于任意复杂的速度介质模型，在地震层析成像等领域得到了广泛应用.但是LTI原算法[[11]存在两个问题：在向前计算节点最小走时时，不能正确追踪逆向传播的射线，相关节点不能得到正确的最小走时[[16-22]；在向后追踪接收点射线路径时，存在不能正确追踪接收点射线路径的可能，算法的稳定性存在不足.文献[17-18]将波前扩展方式与LTI算法基本方程相结合，提出了动态网络最短路径射线追踪算法，该算法在向前计算各节点的最小走时时，从震源点开始，采用波前扩展的方式逐点计算各个节点上的最小走时，改变了LTI原算法的计算方式，确保各节点能得到最小走时；在向后追踪接收点射线路径时，基于互换原理考虑了接收点所有可能的射线路径，确保算法能正确追踪接收点的射线路径.动态网络最短路径射线追踪算法虽然能够解决LTI原算法存在的两个问题，但是其计算效率偏低.文献[22]基于波的传播规律提出了动态网络最短路径射线追踪改进算法，改进并提高了动态网络最短路径射线追踪算法向前计算节点最小走时这一步骤的计算效率，但是，该改进方法的向后追踪方法仍然采用动态网络最短路径射线追踪算法中的方法，计算效率依然偏低.针对动态网络最短路径射线追踪算法中的向后追踪方法存在计算效率低的问题，本文提出了改进方法.首先，在向前计算节点最小走时这一步骤中，不仅计算各节点的最小走时，而且还记录各节点次级源的位置；然后，在向后追踪接收点的射线路径时，利用各节点次级源的位置信息以及波的传播规律对算法进行改进，降低算法的计算量，提高算法的计算效率.1 LTI算法基本方程的推导及接收点文献[8]给出了最短路径射线追踪算法中次级源的确定方法，考虑到LTI算法中的射线可通过单元边界上任意点进行传播，这与最短路径射线追踪算法中射线只能通过节点进行传播不同，因此最短路径射线追踪算法中次级源的确定方法不完全适合LTI算法.为了更好地描述本文提出的改进方法，将LTI算法中次级源的确定方法规定如下：通过接收点所在单元各节段（不包括接收点所在节段）可得到多个接收点走时，其中走时最小值对应的点即为接收点的次级源，如果走时最小值对应多个点，取距离接收点最近的点作为接收点的次级源.如图1所示，若接收点C通过AB节段D点得到的走时小于通过单元其它节段得到的走时，那么D点即为接收点C的次级源. 2 LTI原始算法向后追踪过程存在的问题及动态网络最短路径射线追踪算法存在的不足及改进LTI原始算法在向后追踪接收点射线路径时，首先逐点计算接收点所在单元中各节点走时与节点至接收点的走时之和，然后选出最小的走时之和以及相应的节点，最后将单元中包含该节点的节段作为接收点次级源的可能区域[[11].事实上，接收点的次级源并不一定在这些节段中.此外，LTI原算法在确定接收点次级源的可能区域时，并未排除接收点所在的节段，算法可能会陷入无限循环.因为通过接收点所在节段线性插值得到的接收点走时，可能比通过单元其它节段得到的走时更小，那么接收点取最小走时对应的点可能为接收点本身，算法可能会进入无限循环.以图2所示模型为例，模型尺寸为3 m×3 m，单元大小为1 m×1 m，单元边界划分为2个节段，模型上层、中层以及下层单元的速度分别为525 m/s，515 m/s和505 m/s，震源S的x，y坐标分别为1.5和3.0，接收点R的x，y坐标分别为2.05和0；图中虚线为采用LTI原始算法中的向后追踪方法得到的射线路径，需要说明的是，计算这条射线路径时，已经排除了接收点的次级源位于接收点所在节段的情况，实线为根据动态网络最短路径射线追踪算法中的向后追踪方法[[17-18]得到的射线路径.从图2可以看出，在III号单元中，两种方法的计算结果一样，接收点R的次级源均为R1.在II号单元中，对于LTI原始算法的向后追踪方法，首先须确定II号单元中节点走时与节点至接收点R1走时之和最小的节点，经计算确定为d节点，然后确定cd或de节段为接收点R1次级源的可能区域；对于动态网络最短路径射线追踪算法中的向后追踪方法，接收点R1次级源所在区域为ef节段.计算结果表明：通过ef 节段计算得到的R1点的走时（T=0.005 416 05 s）要小于通过cd或de节段的走时（T=0.005 422 57 s），接收点R1的次级源不在cd或de节段内.造成这一问题的根源在于：LTI算法假定射线路径可以经过单元边界的任意一点，在向后追踪过程中必须考虑所有可能的射线路径，并根据费马原理选择走时最短的那条路径[[11]，而LTI原始算法的具体追踪方法却没有考虑接收点所在单元所有节段的射线，而是采用了一种“简化”方法确定接收点射线路径或次级源的可能区域，由图2所示模型可以看出，这种“简化”方法存在不足.此外，在确定接收点次级源时，若不排除接收点所在的节段，那么在计算新接收点R1的次级源时，由于通过节段cd插值得到的走时，比通过ed节段得到的走时小，计算得到的“次级源”为R1本身，程序将陷入无限循环.动态网络最短路径射线追踪算法中基于互换原理提出的向后追踪方法[[17-18]，考虑了来自接收点所在单元中除接收点所在节段外所有节段的射线，确保了接收点能得到其次级源.但是，该方法计算量大，计算效率低.为了解决这一问题，本文首先在向前计算节点最小走时的步骤中，建立了一个数组，专门用于记录节点次级源的位置信息.然后利用各节点次级源的位置信息以及波的传播规律对动态网络最短路径射线追踪算法中的向后追踪方法进行改进.现以图3所示模型为例，对该向后追踪方法的基本步骤及计算策略进行说明，同时对向后追踪改进方法进行阐述，改进前后的向后追踪方法及计算策略分别如图3和图4所示.具体步骤和分析如下：1）首先将接收点分为以下3种情况，然后对不同的情况采用不同的策略求解接收点的次级源：①接收点位于单元内部，如图3（a）及图4（a）所示的R点.此时利用LTI算法的相关公式计算接收点所在单元各节段至接收点的走时，然后从中选出走时最小值对应的点，如图3（b）及图4（b）中的R1点，这个点就是接收点的次级源.②接收点位于单元边界上，但是非单元边界上的节点.如图3（b），（d）中的R1和R3点，此时除接收点所在的节段外，单元中的其它节段均要计算从该节段至接收点的最小走时，然后从中选出走时最小值对应的点，如图3（c），（e）中的R2和S点.事实上，对于接收点位于单元边界但非单元节点的情况，除非该接收点为原始接收点（图3和图4中的R），否则并不需要对接收点进行全方位的计算.可以采取如下两个步骤确定接收点次级源的可能区域.首先，利用接收点所在节段端点的次级源位置信息确定两个定位点.定位点的确定方法为：若节点端点的次级源不在单元内，或者次级源在单元内但与节点端点处于同一边界，则定位点为节段端点本身，其它情况定位点为节点端点的次级源.如图4（b）所示，接收点R1所在节段的端点为A1和A2，易知接收点R1的次级源位于单元I中，现以端点A1的次级源位置信息为例说明定位点的确定方法，若A1的次级源不在单元I中，比如A1的次级源为单元II中的b点，则由A1的次级源位置信息确定的定位点为A1本身；若A1的次级源与A1处于同一边界，比如次级源位于A1A4或者A1A6边界，则定位点为A1本身；对于其它情况，比如A1的次级源为单元I中的c点，则定位点为A1的次级源c.在图4（b）中，假定A1节点的次级源为A3节点，A2节点的次级源为a1点，则由A1，A2节点的次级源位置信息可以确定两个定位点，分别为A1节点本身以及点a1.然后沿着单元边界连接两个定位点，其中不包含接收点的那条路径即为接收点次级源的可能区域.如图4（b）所示，沿单元边界连接定位点A1和a1可以得到两条路径：A1A3A4a1和A1A6A7a1，其中路径A1A3A4a1不包含接收点R1，因此确定该路径为接收点次级源的可能区域.得到次级源的可能区域后，再确定这个区域内的节段，图4（b）中可能区域A1A3A4a1内的节段为A1A3，A3A4，A4A5，计算从这些节段至接收点的最小走时，其中走时最小的点即为接收点的次级源，如图4（c）中的R2. 同理得到图4（d）中R3点的次级源S.③接收点位于单元边界上，且为单元边界上的节点.如图3（c）中的R2节点，此时需要在R2节点所处的几个单元中执行情况②的计算，如图3（c）.显然，由于在向前计算节点最小走时的过程中，已经记录了各节点的次级源，因此，在改进方法中，可以直接得到接收点的次级源，如图4（c）.2）以步骤1）中获得的次级源为新的接收点，重复步骤1），直至新的接收点为震源.3）依次连接各接收点得到R点的初至波射线路径，如图3（e）和图4（e）所示.对比图3与图4可以看出，改进后的向后追踪方法计算量明显减少.3 数值算例为了对比改进前后向后追踪方法的计算效率，建立了尺寸为2 500 m×600 m的二维模型，如图5所示，其中x值在500～2 000 m之间且y值在200～400 m之间的区域为低速区，速度为500 m/s，其余速度均为4 000 m/s，震源S位于模型上表面的正中间（1 250，0），单元尺寸为5 m×5 m，单元边界划分段数为4，10，20等3种情况，模型下表面的每个单元布置一个接收点，共500个，部分接收点的射线路径如图6所示，分别记录两种向后追踪方法追踪所有接收点的射线路径所耗费的总时间.此外，为了说明向后追踪方法计算效率的提高对整个算法的影响，将改进后的向后追踪方法应用于动态网络最短路径射线追踪改进算法[[22]，然后比较应用改进方法前后动态网络最短路径射线追踪改进算法的总耗时（不包括算法的前处理过程）.计算机CPU主频为3.4 GHz，计算结果如表1所示.由表1可知，向后追踪改进方法的计算效率较高，是改进前向后追踪方法的几倍至几十倍.并且能将动态网络最短路径射线追踪改进算法的计算效率提高1倍左右.为验证向后追踪改进方法对复杂模型的有效性，采用向后追踪改进方法对Marmousi速度模型进行射线追踪，其中节点最小走时的计算采用的是动态网络最短路径射线追踪改进算法.模型尺寸为9 192 m×2 904 m，震源S位于模型左上角（0，0），单元尺寸为24 m×24 m，单元边界划分为10段，共设置5个接收点R1～R5，分别位于模型上表面的3 600，4 800 m，6 000 m，7 200 m和8 400 m，射线追踪结果如图7（a）所示，作为对照，本文给出了单元边界划分为30段时最短路径法的射线追踪结果，如图7（b）所示.两种算法下，各接收点根据射线追踪结果计算的走时如表2所示.计算结果表明，向后追踪改进方法对于复杂的速度模型同样有效.4 结论针对LTI原算法中向后追踪方法存在的无限循环以及可能不能得到正确射线路径的问题，动态网络最短路径射线追踪算法中基于互换原理提出的向后追踪方法能够予以有效的解决.但是该算法存在较多的无效计算.本文根据模型中节点次级源的位置信息以及波的传播规律，提出了改进的向后追踪方法.数值结果表明，改进后的向后追踪方法，其计算效率较之改进前的方法有较大程度的提高；此外，本文提出的向后追踪改进方法能将动态网络最短路径射线追踪改进算法的计算效率提高1倍左右.参考文献[1] JULIAN B R，GUBBINS D.Three-dimensional seismic ray tracing[J].J Geophys，1977，43（1/2）：95-113.[2] 田?h，陈晓非.水平层状介质中的快速两点间射线追踪方法[J]. 地震学报，2005，27（2）：147-154.[3] XU Tao，ZHANG Zhong-jie，GAO Ergen，et al. Segmentally iterative ray tracing in complex 2D and 3D heteroge-neous block models[J]. Bulletin of the Seismological Society of America，2010，100（2）：841-850.[4] VIDALE J.Finite-difference calculation of traveltimes[J].Bulletin of the Seismological Society of America，1988，78（6）：2062-2076.[5] QIN Fu-hao，LUO Yi，OLSENl K B，et al.Finite-difference solution of the eikonal equation along expanding wavefronts[J].Geophysics，1992，57（3）：478-487.[6] 李振春，刘玉莲，张建磊，等.基于矩形网格的有限差分走时计算方法[J]. 地震学报，2004，26（6）：644-650.[7] MOSER T J. Shortest path calculation of seismic rays[J]. Geophysics，1991，56（1）：59-67.[8] 刘洪，孟凡林，李幼铭.计算最小走时和射线路径的界面网全局方法[J].地球物理学报，1995，38（6）：823-832.[9] 赵爱华，徐涛.提高规则网格最短路径方法反射波走时计算精度的走时校正技术[J]. 地球物理学进展，2012，27（5）：1854-1862.[10]BAI Chao-yin，HUANG Guo-jiao，ZHAO Rui. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications[J]. Geophysical Journal International，2010，183（3）：1596-1612. [11]ASAKAWA E，KAWANAKA T. Seismic ray tracing using linear traveltime interpolation[J]. Geophysical Prospecting，1993，41（1）：99-111.[12]赵改善，郝守玲，杨尔皓，等.基于旅行时线性插值的地震射线追踪算法[J].石油物探，1998，37（2）：14-24.[13]CARDARELLI E，CERRETO A. Ray tracing in elliptical anisotropic media using the linear traveltime interpolation （LTI）method applied to traveltime seismic tomography[J].Geophysical Prospecting，2002，50（1）：55-72.[14]聂建新，杨慧珠.地震波旅行时二次/线性联合插值法[J].清华大学学报：自然科学版，2003，43（11）：1495-1498.NIE Jian-xin，YANG Hui-zhu.Quadratic/linear travel timeinterpolationof seismic ray tracing[J].J Tsinghua Univ ：Sci&Tech，2003，43（11）：1495-1498.（In Chinese）[15]ZHANG Jian-zhong，HUANG Yue-qin，SONG Lin-ping，et al.Fast and accurate 3-D ray tracing using bilinear traveltime interpolation and the wave front group marching[J]. Geophysical Journal International，2011，184（3）：1327-1340.[16]黄靓，黄政宇.线性插值射线追踪的改进方法[J]. 湘潭大学自然科学学报，2002，24（4）：105-108.[17]张建中，陈世军，余大祥.最短路径射线追踪方法及其改进[J].地球物理学进展，2003，18（1）：146-150.[18]张建中，陈世军，徐初伟.动态网络最短路径射线追踪[J]. 地球物理学报，2004，47（5）：899-904.[19]黄靓.混凝土超声波层析成像的理论方法和试验研究[D]. 长沙：湖南大学土木工程学院，2008：33-36.[20]张东，谢宝莲，杨艳，等.一种改进的线性走时插值射线追踪算法[J]. 地球物理学报，2009，52（1）：200-205.[21]卢江波，方志.线性走时插值射线追踪算法的改进[J].湖南大学学报：自然科学版，2014，41（1）：39-44.[22]卢江波，方志.一种线性走时插值射线追踪改进算法[J].地震学报，2014，36（6）：1089-1100.。

最长公共子序列问题的改进快速算法
李欣;舒风笛
【期刊名称】《计算机应用研究》
【年(卷),期】2000(017)002
【摘要】现在几个最常用的解决最长公共子序列(LCS)问题的算法的时间复杂度分别是O(pn),O(n(m-p)).这里m、n为两个待比较字符串的长度,p是最长公共子串的长度.给出一种时间复杂度为O(p(m-p)),空间复杂度为O(m+n)的算法.与以前的算法相比,不管在p<<m的情况下,还是在p接近m时,这种算法都有更快的速度.【总页数】3页(P28-30)
【作者】李欣;舒风笛
【作者单位】北京大学软件工程实验室,北京,100871;武汉大学计科系,武
汉,430072
【正文语种】中文
【中图分类】TP3
【相关文献】
1.并行改进回溯算法实现N皇后问题的快速计数 [J], 韩宇南;吕英华;黄小红
2.流水作业调度问题的快速进入启发式算法改进 [J], 王崐
3.流水作业调度问题的快速进入启发式算法改进 [J], 王崐;;
4.基于改进蚁群算法快速求解迷宫路径问题研究 [J], 杨乐;向凤红;毛剑琳
5.基于改进蚁群算法快速求解迷宫路径问题研究 [J], 杨乐;向凤红;毛剑琳
因版权原因，仅展示原文概要，查看原文内容请购买。