使用Playground学习数值算法

本文由CocoaChina翻译组成员leon（社区ID）翻译自raywenderlich
原文：Numerical Algorithms using Playgrounds

中学时，没有什么比数学图纸更恐怖的了。

许多问题没有现成解析方案，还有一些问题虽然可以解决，但需要大量的计算工作。这种情况下，你需要算法来辅助。
希望你没被算法整吐。万一你真的吐了，也要看到好的一面：又可以再吃顿午餐了:]
在这个教程中，你将学习到算法的概念，以及如何使用它们来解决难以分析的问题。通过playground，很容易使解决方案可视化，易于查看。
即使你不是个数学爱好者，对物理或计算机科学也不大感兴趣，你仍会在这个教程中找到一些有意思的东西。你只需要有一些微积分和基础物理学的知识。
你将学到如何使用数值算法解决两个困难问题。但是为了讲的更清楚点，这两个问题也可通过解析法解决。算法在解析法无法工作时更为理想，即便如此，通过比较两种方法，也能更加容易的理解它们的本质。
数值算法是什么？
简单说来，数值算法是一类解决数学问题的方法，它们不依赖于闭式解析解。
问题来了，什么是闭式解析解？
若有一个公式，我们可以使用已知数，通过有限的操作步骤，可以获得一个准确的结果，即称之为闭式解析解。
再简单点来理解一下：如果你可以使用代数的方式，找到一个表达式来解决一个未知量问题，代入某些已知数即可得到结果，就说明你已经得到了一个闭式解析方法。
何时使用数值算法？
许多问题没有现成解析方案，还有一些问题虽然可以解决，但需要大量的计算工作。这种情况下，你需要算法来辅助。
例如，假定你需要编写一个物理引擎，用来计算许多物体在有限时间内的行为。这时你就可以使用数值算法来更快地得到结果。
这样做也有副作用，更快的计算速度意味着结果不会十分精确。然而，大多数情况下，这样的结果已经够用了。
天气预报就得益于数值计算。天气变化的速度、影响因素的数量都是十分惊人的。这是一个非常复杂的系统，只有数值模拟才能完成预知未来的任务。
可能正是因为缺乏这些算法，你的iPhone总是告诉你要下雨了，而外面还是艳阳高照。


开始
作为热身活动，我们来玩个游戏，接下来，你将计算出一个给定数字的平方根。这两个方法都将使用二分法（bisection method）。神奇的是，你可能已经知道了这个方法，但是不知道它的名字。
回想一下，在儿童时代，我们可能玩过这样的游戏：选中100以内的一个数字，另外一个人来猜。你只会提示他猜的数字大了或者小了。
游戏开始。小明告诉小强开始猜，小强说我猜1，小明说小了，小强又猜2，小明说小了。小强接下来再猜3,4...最后选中了5，终于对了。
5步就猜中了，不错。但是如果小明选的是78，猜中就需要花点时间。
这个游戏如果使用二分法（bisection method）来玩的话，猜中的速度会快很多。
二分法
我们知道数字介于1和100之间，因此除了一个一个的猜，或者随便瞎猜外。我们把数字分为两个区间：a = [1,50], b = [51, 100]。
接下来我们判断数字是介于哪一个区间，这很简单，只用问数字是不是50。如果数字小于50，那么区间b就不用考虑了。接下来我们继续把a区间分成两半，再重复这个步骤。
例如：
数字是60.区间为：a = [1,50], b = [51, 100]。
第一步，小强说50（也即是a区间的上限），小明说小了。这时小强就知道了数字肯定在b区间上。
第二步，分解b区间为： c=[51,75]，d=[76,100]。还是选择c区间的上限75，小强告诉他大了。这就意味着数字肯定在c区间上。
继续分解。。。
使用这个方法，7次尝试即可得到正确答案，一个一个试则需要60次。
1. 50 -> 小了
2. 75 -> 大了
3. 62 -> 大了
4. 56 -> 小了
5. 59 -> 小了
6. 61 -> 大了
7. 60 -> 对了！
计算x的平方根，过程也类似。数字x的平方根介于0和x之间。也可以记做：(0,x]。如果数字x大于或等于1，可以记做：[1, x]。
分解区间，得到a=(0, x/2]，b=(x/2, x]。
如果数字x是9，区间是[1, 9]，分解后的区间为a=[1, 5]，b=(5, 9]，中间值m为(1+9)/2 = 5。
下一步，检查m * m - x，是否大于期望的精度。在这里，我们检查m * m 大于或小于等于x。如果大于，我们使用(0, x/2]区间继续检查，否则，使用(x/2, x]区间。
看一下实际步骤，初始化m=5, 期望准确值为0.1:
1. 计算m * m - x: 5 * 5 - 9 = 11
2. 检查结果是否小于等于期望值：25 - 11 <= 0.1?
3. 显然不满足，继续下一个区间：m * m是否大于9？是。接下来使用区间[1, 5]，测试值为(1+5)/2=3。
4. 计算m * m - x：３* 3 - 9 = 0
5. 查查是否满足期望值：9 - 9 <= 0.1?
6. 搞定。
备注：你可能想知道括号是什么意思？区间有固定的格式（下限和上限）。不同的括号代表不同的区间边界。圆括号表示边界不在区间范围内（即开区间），而方括号表示边界在区间范围内（闭区间）。如(0, a] 包含a而不包含0.在上面的例子中：区间 (0, a/2]包含a/2而不包含0；区间(a/2, a]表示所有大于a/2, 小于a的数。
在Playground上使用二分法
现在，是时候应用学到的理论了。我们来自己实现二分算法。创建一个新的playground文件，添加如下的代码：

func bisection(x: Double) -> Double {
//1
    var lower = 1.0
    //2
    var upper = x
    //3
    var m = (lower + upper) / 2
    var epsilon = 1e-10
    //4
    while (fabs(m * m - x) > epsilon) {
    //5
        m = (lower + upper) / 2
        if m * m > x {
            upper = m
        } else {
            lower = m
        }
    }
 
    return m
}

看一下代码各部分的含义：
1. 设置区间下限为1
2. 设置区间上限为x
3. 找到中间值，并定义期望精确度
4. 检查操作是否能满足精确度
5. 如果不满足，找到新的中间值，定义新的区间，继续查找。
添加如下代码来测试该函数：

let bis = bisection(2.5)

在 m = (lower + upper) / 2这一行的右边，可以看到代码执行了35次，意味着找到结果需要35步。
马上我们就能看到playground一个可爱功能带来的好处：数值的完整历史都可以查看。
二分法可以成功的算出实际结果的近似值。通过数值的历史记录图，就可以看到数值算法是如何逼近正确的解。
按下option+cmd+enter打开辅助编辑器，点击m = (lower + upper) / 2代码行右边的圆按钮在辅助编辑器中添加历史记录。


你会看到方法一点点的跳转到正确结果上。

古典数学也搞的定
下一个算法需要追溯到古代，它起源于公元前1750年的古巴比伦，在公元100年前亚历山大的希罗所著《度量论》（Metrica ）一书中有所描述。这就是为何它被称作“希罗方法”。可以通过这个 链接 了解更多关于希罗的知识。
这个方法使用函数



，这里你想要计算a的平方根。如果你能找到函数曲线的“0值点（zero point）”，这个点上的函数值为0，那么求出x值就能得到a的平方根。
要完成这项工作，我们首先选择任意x值作为起始值，计算该值对应点的tangent值（函数的切线），然后查看tangent线上的0点（即该直线和x轴的交点）。然后我们使用这个0点再次作为起始值，重复多次直到满足精度。
因为每计算一次tangent值，都会更加逼近真正的0值，
这个过程会逐渐逼近真正的结果。下图展示了求解a=9时，a的平方根，且起始值为1.
起始点x0 = 1，生成一条红色的tangent线，接着生成了x1点，又生成了紫色的线；又生成了x2点，连接了蓝色的线，最终找到了正确答案。


当古典数学邂逅了Playground
我们有很多古巴比伦人没有的好东西，比如：playground。在playground上添加如下代码，并查看：

func heron(x: Double) -> Double {
  //1
  var xOld = 0.0
  var xNew = (x + 1.0) / 2.0
  var epsilon = 1e-10
 
  //2
  while (fabs(xNew - xOld) > epsilon) {
    //3
    xOld = xNew
    xNew = (xOld + x / xOld) / 2
  }
 
  return xNew
}

这段代码做了什么？
1. 创建一个变量来存储当前结果。xOld是上次计算的结果，而xNew是实际结果。
找到初始化xNew的方法是一个良好的起点

epsilon是期望的精确度

这个例子中，我们计算平方根精确度为小数点后10位。

2. 使用while循环查找是否达到期望的精确度。
3. 如果达不到精确度，设置xNew为xOld，继续下一次迭代。
使用如下代码验证该函数的作用：

let her = heron(2.5)

希罗方法只需要5次迭代即可找到正确结果。
在代码行xNew = (xOld + x/xOld)/2右边点击圆角button，添加值历史，就能看到第一次迭代就能找到一个不错的近似值。


模拟谐振子运动

在这个章节中，我们来看看如何使用数值积分算法来模拟简谐系统运动——一种基本的动态系统。
这个系统可以描述很多现象，比如钟摆的摆动、弹簧的振动。特别说来，它可以描述某段时间内物体移动了一定偏移量的场景。
这个例子中，假定有一个有质量的物体连接在弹簧上。为了简化问题，我们忽略掉阻力和重力。因此唯一的作用力就是弹簧的弹力，它将物体拉回到原来的位置。
在这样的假定下，你只需要使用两个力学定律来处理问题：
牛顿第二运动定律

, 它描述了物体上的作用力和加速度的关系。

胡克定律

，它描述了弹力和物体偏移量之间的比例关系。这里k是弹力系数而x是物体的偏移量。

简振方程
因为弹力是物体上唯一的作用力，我们将上面的两个方程式整理：



简化后写作：




也被记作

，也即是共振频率的平方。
更精确的方程式如下：



其中A是振幅，这里即是物体的偏移量，

被称为相位差。这两个值初始化时都被设置为定值。
如果设置时间t=0，则

，且

，你可以简单的计算出振幅和相位差：





来看一个例子，我们有一个质量为2kg的物体连接在弹簧上，弹簧的弹力系数为k=196N/m。在初始时间t=0时，弹簧移动了0.1米。我们可以通过公式计算振幅

、相差

和共振频率

。







在Playground上实验一下：
使用该公式计算任意时间点对应的值之前，我们需要添加一些代码。
回到playground文件，在最后添加如下代码：

//1
typealias Solver = (Double, Double, Double, Double, Double) -> Void
 
//2
struct HarmonicOscillator {
  var kSpring = 0.0
  var mass = 0.0
  var phase = 0.0
  var amplitude = 0.0
  var deltaT = 0.0
 
  init(kSpring: Double, mass: Double, phase: Double, amplitude: Double, deltaT: Double) {
    self.kSpring = kSpring
    self.mass = mass
    self.phase = phase
    self.amplitude = amplitude
    self.deltaT = deltaT
  }
 
  //3
  func solveUsingSolver(solver: Solver) {
    solver(kSpring, mass, phase, amplitude, deltaT)
  }
}

这些代码块中做了什么？
1. 定义一个函数类型的的typealias，函数有5个Double类型的参数，返回为空。
2. 创建一个struct来定义谐振子。
3. 这个函数只是简单的创建用来解振动问题的Clousure。（并无真实的计算代码）
再来看一下精确方案
代码如下：

func solveExact(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) {
  var x = 0.0
  //1
  let omega = sqrt(kSpring / mass)
 
  var i = 0.0
 
  while i < 100.0 {
//2
    x =  amplitude * sin(omega * i + phase)
    i += t
  }
}

这个方法包含了所有解决运动方程需要的参数。
1. 计算共振频率
2. 在while循环中计算物体当前的位置，i用作下次计算的自增变量。
添加如下测试代码：

let osci = HarmonicOscillator(kSpring: 0.5, mass: 10, phase: 10, amplitude: 50, deltaT: 0.1)
osci.solveUsingSolver(solveExact)

这个方案中的方法有点奇怪：它有参数传入，但不返回数据，也没有显示任何东西。
这样做有什么好处？
使用这个函数的目的在于，在while循环中，可以算出振动过程中具体的动态值。在Playground中，可以观察这些动态值的历史记录。
在x = amplitude * sin(omega * i + phase) 处添加值历史记录，我们就能看到运动的轨迹。


既然第一个精确算法已经验证通过，下面我们看一下数值计算方案。

欧拉方法
欧拉方法是用来求数值积分的一种方法。该方法1768年再Leonhard Euler所著Institutiones Calculi Integralis《积分学原理》一书中首次提出。要查看该方法的详情，请参考：http://en.wikipedia.org/wiki/Euler_method
欧拉方法的核心思想是通过使用折线逼近曲线。
具体方法是计算一个给定点的斜率，然后绘制一条同样斜率的折线。在这条折线的末尾，继续计算斜率，绘制另外一条线。正如你所见，该算法的精确度取决于折线的长度。
你想知道deltaT做了什么吗？
该数值算法中需要使用一个步长，这对算法的精确度很重要。大步长导致低精确度，但是执行速度块。反之，精确度会提高，速度会降低。
deltaT变量就代表了步长(step size)。我们初始化这个值为0.1， 意味着我们计算每0.1秒物体移动的位置。在欧拉算法中，意味着折线在x轴上的投影长度为0.1。


在编写代码前，你需要再看一下这个公式：



二阶微分方程可以转化为两个一阶微分方程。

可以被写为

和


我们可以用差商来完成转换：



也会得到：



有了这些等式，我们可以直接实现欧拉方法。
在solveExact方法后添加代码：

func solveEuler(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) {
  //1
  var x = amplitude * sin(phase)
  let omega = sqrt(kSpring / mass)
  var i = 0.0
  //2
  var v = amplitude * omega * cos(phase);
  var vold = v
  var xoldEuler = x
 
  while i < 100 {
    //3
    v -= omega * omega  * x * t
    //4
    x += vold * t
    xoldEuler = x
 
    vold = v
    i+=t
    }
}

代码的内容：
1. 设置当前的位置，和omega的值。
2. 计算当前速度。
3. 在while循环中，通过一阶函数计算出新的速度。
4. 使用速度计算新的位置，在结束处，使用步长t增加i的值。
现在，添加以下代码测试：

osci.solveUsingSolver(solveEuler)

在xoldEuler = x位置添加值历史，并查看。你会看到这个方法显示一个振幅增加的正弦函数。因为欧拉方法并不是一个精确算法，而这里过大的步长0.1也导致了计算结果不精确。


以下是另外一个函数图像，步长为0.01，这样明显更好。因此，使用欧拉方法时你要想到，步长越小，结果越好。但是使用折中的步长，执行起来更为简单。


速度Verlet算法(Velocity Verlet)
最后讨论的算法叫做速度Verlet。它和欧拉方法的思路一样，但是计算新位置的方式有些许差异。
欧拉在计算新位置时，忽略实际的加速度，公式为：

,而速度Verlet算法在计算时考虑到了加速度：

。这再步长相等的时候，结果更优。
在solveEuler方法后添加新的代码：

func solveVerlet(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) {
  //1
  var x = amplitude * sin(phase)
  var xoldVerlet = x
  let omega = sqrt(kSpring / mass)
 
  var v = amplitude * omega * cos(phase)
  var vold = v
  var i = 0.0
 
  while i < 100 {
    //2
    x = xoldVerlet + v * t + 0.5 * omega * omega * t * t
    v -= omega * omega * x * t
    xoldVerlet = x
    i+=t
  }
}

代码的内容：
1. 循环前的所有代码和欧拉方法中的一样。
2. 根据速度计算出新的位置，增加i的值。
在Playground中测试代码:

osci.solveUsingSolver(solveVerlet)

别忘了在xoldVerlet = x代码行上添加值历史记录：


接下来做什么？
你可以通过 此链接 下载完整工程.
希望你在这场数值世界旅行中获得乐趣。了解一些算法，即便只是一些古典数学的趣味历史，都会对你带来帮助。
继续探索之前，你可以再改改deltaT，看看有什么不同的结果。但是要记住playground速度并不快，选择一个小的步长可能导致一段时间内Xcode无法操作。在我的Macbook上，设置deltaT到0.01就让我的机器定住了好几分钟。
使用你学到的算法，可以解决许多问题，比如说：n-body 问题。
如果你觉得自己的算法基础不牢固。可以在YouTube上看看MIT的 算法简介课程。
最后，有时间可以在本站点查看更多关于这个课题的教程。
如果你对于使用数值算法有任何的疑问或建议，可以在下面评论。期望看到你的想法，并讨论在app中可以用到的很酷的数学知识。