使用Playground学习数值算法

移动开发 算法
在这个教程中,你将学习到算法的概念,以及如何使用它们来解决难以分析的问题。通过playground,很容易使解决方案可视化,易于查看。

chart-250x250.jpg

 

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

许多问题没有现成解析方案,还有一些问题虽然可以解决,但需要大量的计算工作。这种情况下,你需要算法来辅助。

希望你没被算法整吐。万一你真的吐了,也要看到好的一面:又可以再吃顿午餐了:]

在这个教程中,你将学习到算法的概念,以及如何使用它们来解决难以分析的问题。通过playground,很容易使解决方案可视化,易于查看。

即使你不是个数学爱好者,对物理或计算机科学也不大感兴趣,你仍会在这个教程中找到一些有意思的东西。你只需要有一些微积分和基础物理学的知识。

你将学到如何使用数值算法解决两个困难问题。但是为了讲的更清楚点,这两个问题也可通过解析法解决。算法在解析法无法工作时更为理想,即便如此,通过比较两种方法,也能更加容易的理解它们的本质。

数值算法是什么?

简单说来,数值算法是一类解决数学问题的方法,它们不依赖于闭式解析解。

问题来了,什么是闭式解析解?

若有一个公式,我们可以使用已知数,通过有限的操作步骤,可以获得一个准确的结果,即称之为闭式解析解。

再简单点来理解一下:如果你可以使用代数的方式,找到一个表达式来解决一个未知量问题,代入某些已知数即可得到结果,就说明你已经得到了一个闭式解析方法。

何时使用数值算法?

许多问题没有现成解析方案,还有一些问题虽然可以解决,但需要大量的计算工作。这种情况下,你需要算法来辅助。

例如,假定你需要编写一个物理引擎,用来计算许多物体在有限时间内的行为。这时你就可以使用数值算法来更快地得到结果。

这样做也有副作用,更快的计算速度意味着结果不会十分精确。然而,大多数情况下,这样的结果已经够用了。

天气预报就得益于数值计算。天气变化的速度、影响因素的数量都是十分惊人的。这是一个非常复杂的系统,只有数值模拟才能完成预知未来的任务。

可能正是因为缺乏这些算法,你的iPhone总是告诉你要下雨了,而外面还是艳阳高照。

[[142780]]

开始

作为热身活动,我们来玩个游戏,接下来,你将计算出一个给定数字的平方根。这两个方法都将使用二分法(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* 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文件,添加如下的代码:

  1. func bisection(x: Double) -> Double { 
  2. //1 
  3.     var lower = 1.0 
  4.     //2 
  5.     var upper = x 
  6.     //3 
  7.     var m = (lower + upper) / 2 
  8.     var epsilon = 1e-10 
  9.     //4 
  10.     while (fabs(m * m - x) > epsilon) { 
  11.     //5 
  12.         m = (lower + upper) / 2 
  13.         if m * m > x { 
  14.             upper = m 
  15.         } else { 
  16.             lower = m 
  17.         } 
  18.     } 
  19.    
  20.     return m 
看一下代码各部分的含义:

1. 设置区间下限为1

2. 设置区间上限为x

3. 找到中间值,并定义期望精确度

4. 检查操作是否能满足精确度

5. 如果不满足,找到新的中间值,定义新的区间,继续查找。

添加如下代码来测试该函数:

  1. let bis = bisection(2.5

在 m = (lower + upper) / 2这一行的右边,可以看到代码执行了35次,意味着找到结果需要35步。

#p#

马上我们就能看到playground一个可爱功能带来的好处:数值的完整历史都可以查看。

二分法可以成功的算出实际结果的近似值。通过数值的历史记录图,就可以看到数值算法是如何逼近正确的解。

按下option+cmd+enter打开辅助编辑器,点击m = (lower + upper) / 2代码行右边的圆按钮在辅助编辑器中添加历史记录。

Screen-Shot-2015-05-13-at-11.07.08-AM1-603x500.png

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

古典数学也搞的定

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

这个方法使用函数latex.pnglatex.png,这里你想要计算a的平方根。如果你能找到函数曲线的“0值点(zero point)”,这个点上的函数值为0,那么求出x值就能得到a的平方根。

要完成这项工作,我们首先选择任意x值作为起始值,计算该值对应点的tangent值(函数的切线),然后查看tangent线上的0点(即该直线和x轴的交点)。然后我们使用这个0点再次作为起始值,重复多次直到满足精度。

因为每计算一次tangent值,都会更加逼近真正的0值,

这个过程会逐渐逼近真正的结果。下图展示了求解a=9时,a的平方根,且起始值为1.

起始点x0 = 1,生成一条红色的tangent线,接着生成了x1点,又生成了紫色的线;又生成了x2点,连接了蓝色的线,最终找到了正确答案。

Heron.jpg

当古典数学邂逅了Playground

我们有很多古巴比伦人没有的好东西,比如:playground。在playground上添加如下代码,并查看:

  1. func heron(x: Double) -> Double { 
  2.   //1 
  3.   var xOld = 0.0 
  4.   var xNew = (x + 1.0) / 2.0 
  5.   var epsilon = 1e-10 
  6.    
  7.   //2 
  8.   while (fabs(xNew - xOld) > epsilon) { 
  9.     //3 
  10.     xOld = xNew 
  11.     xNew = (xOld + x / xOld) / 2 
  12.   } 
  13.    
  14.   return xNew 

这段代码做了什么?

1. 创建一个变量来存储当前结果。xOld是上次计算的结果,而xNew是实际结果。

  • 找到初始化xNew的方法是一个良好的起点

  • epsilon是期望的精确度

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

2. 使用while循环查找是否达到期望的精确度。

3. 如果达不到精确度,设置xNew为xOld,继续下一次迭代。

使用如下代码验证该函数的作用:

  1.      
  2. let her = heron(2.5

希罗方法只需要5次迭代即可找到正确结果。

在代码行xNew = (xOld + x/xOld)/2右边点击圆角button,添加值历史,就能看到第一次迭代就能找到一个不错的近似值。

Heron.png

模拟谐振子运动

在这个章节中,我们来看看如何使用数值积分算法来模拟简谐系统运动——一种基本的动态系统。

这个系统可以描述很多现象,比如钟摆的摆动、弹簧的振动。特别说来,它可以描述某段时间内物体移动了一定偏移量的场景。

这个例子中,假定有一个有质量的物体连接在弹簧上。为了简化问题,我们忽略掉阻力和重力。因此唯一的作用力就是弹簧的弹力,它将物体拉回到原来的位置。

在这样的假定下,你只需要使用两个力学定律来处理问题:

  • 牛顿第二运动定律latex (1).png, 它描述了物体上的作用力和加速度的关系。

  • 胡克定律latex (2).png,它描述了弹力和物体偏移量之间的比例关系。这里k是弹力系数而x是物体的偏移量。

简振方程

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

01.png

简化后写作:

02.png

05.png也被记作04.png,也即是共振频率的平方。

更精确的方程式如下:

033.png

其中A是振幅,这里即是物体的偏移量,06.png被称为相位差。这两个值初始化时都被设置为定值。

如果设置时间t=0,则07.png,且08.png,你可以简单的计算出振幅和相位差:

09.png
10.png

来看一个例子,我们有一个质量为2kg的物体连接在弹簧上,弹簧的弹力系数为k=196N/m。在初始时间t=0时,弹簧移动了0.1米。我们可以通过公式计算振幅11.png、相差12.png和共振频率13.png

14.png
15.png
16.png

在Playground上实验一下:

使用该公式计算任意时间点对应的值之前,我们需要添加一些代码。

回到playground文件,在最后添加如下代码:

 
  1. //1 
  2. typealias Solver = (Double, Double, Double, Double, Double) -> Void 
  3.    
  4. //2 
  5. struct HarmonicOscillator { 
  6.   var kSpring = 0.0 
  7.   var mass = 0.0 
  8.   var phase = 0.0 
  9.   var amplitude = 0.0 
  10.   var deltaT = 0.0 
  11.    
  12.   init(kSpring: Double, mass: Double, phase: Double, amplitude: Double, deltaT: Double) { 
  13.     self.kSpring = kSpring 
  14.     self.mass = mass 
  15.     self.phase = phase 
  16.     self.amplitude = amplitude 
  17.     self.deltaT = deltaT 
  18.   } 
  19.    
  20.   //3 
  21.   func solveUsingSolver(solver: Solver) { 
  22.     solver(kSpring, mass, phase, amplitude, deltaT) 
  23.   } 

这些代码块中做了什么?

1. 定义一个函数类型的的typealias,函数有5个Double类型的参数,返回为空。

2. 创建一个struct来定义谐振子。

3. 这个函数只是简单的创建用来解振动问题的Clousure。(并无真实的计算代码)

再来看一下精确方案

代码如下:

 
  1. func solveExact(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) { 
  2.   var x = 0.0 
  3.   //1 
  4.   let omega = sqrt(kSpring / mass) 
  5.    
  6.   var i = 0.0 
  7.    
  8.   while i < 100.0 { 
  9. //2 
  10.     x =  amplitude * sin(omega * i + phase) 
  11.     i += t 
  12.   } 

这个方法包含了所有解决运动方程需要的参数。

1. 计算共振频率

2. 在while循环中计算物体当前的位置,i用作下次计算的自增变量。

添加如下测试代码:

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

这个方案中的方法有点奇怪:它有参数传入,但不返回数据,也没有显示任何东西。

#p#

这样做有什么好处?

使用这个函数的目的在于,在while循环中,可以算出振动过程中具体的动态值。在Playground中,可以观察这些动态值的历史记录。

在x = amplitude * sin(omega * i + phase) 处添加值历史记录,我们就能看到运动的轨迹。

[[142803]]

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

欧拉方法

欧拉方法是用来求数值积分的一种方法。该方法1768年再Leonhard Euler所著Institutiones Calculi Integralis《积分学原理》一书中首次提出。要查看该方法的详情,请参考:http://en.wikipedia.org/wiki/Euler_method

欧拉方法的核心思想是通过使用折线逼近曲线。

具体方法是计算一个给定点的斜率,然后绘制一条同样斜率的折线。在这条折线的末尾,继续计算斜率,绘制另外一条线。正如你所见,该算法的精确度取决于折线的长度。

你想知道deltaT做了什么吗?

该数值算法中需要使用一个步长,这对算法的精确度很重要。大步长导致低精确度,但是执行速度块。反之,精确度会提高,速度会降低。

deltaT变量就代表了步长(step size)。我们初始化这个值为0.1, 意味着我们计算每0.1秒物体移动的位置。在欧拉算法中,意味着折线在x轴上的投影长度为0.1。

Euler_Idea.png

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

17.png

二阶微分方程可以转化为两个一阶微分方程。18.png可以被写为19.png20.png

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

21.png

也会得到:

22.png

有了这些等式,我们可以直接实现欧拉方法。

在solveExact方法后添加代码:

 
  1. func solveEuler(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) { 
  2.   //1 
  3.   var x = amplitude * sin(phase) 
  4.   let omega = sqrt(kSpring / mass) 
  5.   var i = 0.0 
  6.   //2 
  7.   var v = amplitude * omega * cos(phase); 
  8.   var vold = v 
  9.   var xoldEuler = x 
  10.    
  11.   while i < 100 { 
  12.     //3 
  13.     v -= omega * omega  * x * t 
  14.     //4 
  15.     x += vold * t 
  16.     xoldEuler = x 
  17.    
  18.     vold = v 
  19.     i+=t 
  20.     } 

代码的内容:

1. 设置当前的位置,和omega的值。

2. 计算当前速度。

3. 在while循环中,通过一阶函数计算出新的速度。

4. 使用速度计算新的位置,在结束处,使用步长t增加i的值。

现在,添加以下代码测试:

  1. osci.solveUsingSolver(solveEuler) 

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

[[142811]]

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

[[142812]]

速度Verlet算法(Velocity Verlet)

最后讨论的算法叫做速度Verlet。它和欧拉方法的思路一样,但是计算新位置的方式有些许差异。

欧拉在计算新位置时,忽略实际的加速度,公式为:25.png,而速度Verlet算法在计算时考虑到了加速度:26.png。这再步长相等的时候,结果更优。

在solveEuler方法后添加新的代码:

 
  1. func solveVerlet(amplitude: Double, phase: Double, kSpring: Double, mass: Double, t: Double) { 
  2.   //1 
  3.   var x = amplitude * sin(phase) 
  4.   var xoldVerlet = x 
  5.   let omega = sqrt(kSpring / mass) 
  6.    
  7.   var v = amplitude * omega * cos(phase) 
  8.   var vold = v 
  9.   var i = 0.0 
  10.    
  11.   while i < 100 { 
  12.     //2 
  13.     x = xoldVerlet + v * t + 0.5 * omega * omega * t * t 
  14.     v -= omega * omega * x * t 
  15.     xoldVerlet = x 
  16.     i+=t 
  17.   } 

代码的内容:

1. 循环前的所有代码和欧拉方法中的一样。

2. 根据速度计算出新的位置,增加i的值。

在Playground中测试代码:

  1. osci.solveUsingSolver(solveVerlet) 
别忘了在xoldVerlet = x代码行上添加值历史记录:

[[142815]]

接下来做什么?

你可以通过 此链接 下载完整工程.

希望你在这场数值世界旅行中获得乐趣。了解一些算法,即便只是一些古典数学的趣味历史,都会对你带来帮助。

责任编辑:倪明
相关推荐

2022-04-06 10:06:37

判断算法数值校验

2010-01-07 17:58:49

JSON数值

2017-07-14 10:35:06

2021-07-05 06:39:59

经典算法无序数组

2021-04-16 11:31:24

人工智能深度学习

2018-04-28 16:20:31

机器学习算法分发链路

2020-12-16 15:56:26

机器学习人工智能Python

2009-03-11 10:44:49

.netvb.netArray

2020-06-18 16:05:20

机器学习人工智能算法

2016-11-15 15:02:00

机器学习算法

2021-02-01 10:17:14

编程C语言计算机

2023-07-28 08:13:30

2020-02-03 08:00:00

机器学习人工智能AI

2017-07-13 10:03:43

优化算法Adam深度学习

2020-12-11 12:52:58

Java技术开发

2010-06-01 15:16:36

Zabbix使用

2023-05-16 12:51:52

人工智能GPT-3

2019-03-20 07:50:47

机器学习算法线性回归

2018-01-09 13:42:37

集成学习算法

2014-06-17 09:55:24

机器学习
点赞
收藏

51CTO技术栈公众号