重写牛顿法的例子来使用函数
上节课的最后一个例子显然是介绍函数的一个候选例子。这个例子用到了两个函数，一个是多项式f(x)我们要求的是它的根，另一个是它的导数。
下面是用这两个函数重写的例子。
#程序的根f (x) = x ^ 4 + 2 x ^ 3 - 14 x ^ 2 + 2 + 1 #通过牛顿法#定义函数及其导数def f (x):返回(((x + 2) * x-14) * x + 2) * x + 1 def fp (x):返回((4 * x + 6) * x-28) * x + 2 =输入(“输入一个起始猜根:“)x =浮动(x) y = f (x)公差= 10 e-6 y <)或y >公差:#计算f (x) yp = fp (x) #做一个轮牛顿法x = x - y / yp y = f (x)打印(“根估计”+ str (x))
构造一个有用的函数：计算exp(x)
大多数编程系统都包含有用的函数库，Python也不例外。例如，要使用标准的数学函数，如exp(x)或sin(x)，我们可以导入math模块。
>>>导入数学>>>打印（str(math.exp(3.3))） 27.112638920657883
作为编写函数的练习，让我们看看是否可以构造自己的函数来计算exp(x)。这样做的关键是exp(x)的幂级数表达式：
这是一个可以计算这个幂级数的函数的初稿。
def exp1(x): ““”第一次尝试计算exp(x)。”“” sum = 1对于n在范围（1,16）：sum += x**n/math.factorial(n)返回sum
由于幂级数的和是一个无限的和，我们将不得不满足于一个有限的近似。我通过试错法确定了n = 15的截止值，因为这样可以用最少的求和项产生最好的结果。下面是一个简单测试程序的一些输出，它证实了这一点。
Math.exp (0.5)= 2.0137527074704766 exp1(0.5)= 2.0137527074704766 . (0.3)= 1.3498588075760032
接下来，我们介绍一个简单的优化。由于x**n和math.factorial(n)的计算成本都很高，因此我们利用了求和中连续项之间的简单关系。
这意味着我们可以很容易地从它之前的项计算和中的每一个新项。
下面是使用这种优化的exp函数的一个版本。
defexp2 (x): ““”指数函数的更精细的版本。”“” sum = 1 + x term = x，对于range（2,16）中的n: term = term*x/n sum += term返回sum
下面是一个快速的测试运行，以确认这当然产生相同的数值结果与早期的，效率较低的版本。
Math.exp (0.3)= 1.3498588075760032 exp1(0.3)= 1.349858807576003 exp2(0.3)= 1.349858807576003 exp2(0.3)= 2.0137527074704766 exp1(0.5)= 2.0137527074704766 exp2(0.5)= 2.0137527074704766
我们的指数计算函数似乎工作得很好，但它还有最后一个严重的弱点。幂级数近似的问题在于，它们对赢博体育x值的收敛性都不一样。这里有另一个测试表明，当我们代入一个更大的x时，近似会变得更差。
Math.exp (3.3)= 27.112638920657883
当输入x变大时，我们可以通过增加和中的项数来解决这个问题。相反，我要采取另一种策略。由于exp2对于输入x在0到1的范围内似乎做得很好，所以技巧是将给定的任何x映射到该范围内。我们可以通过以下策略做到这一点：
- 计算使x/2n < 1的最小幂n。
- 计算exp (x / 2 n)。
- 请注意,
- 我们可以通过重复平方计算后一个表达式：
Exp (x) = Exp (x/2n))2n
（exp(x/2n)）2n = (exp(x/2n))2)2n-1
现在是指数计算函数的最终版本它利用了这个策略。
def exp3(x): "“ exp的最终版本-适用于更大范围的x。”“” #最初我们使用x的绝对值z = math.fabs(x) #计算最小的n，使z/2**n <= 1 n = 1，而z > 1: n += 1 z/ = 2.0 #因为z现在在0 <= z <= 1的范围内，exp2将做得很好结果= exp2(z) #因为现在结果= exp(z/2**n)，为了计算exp(z)我们现在需要计算# result**(2**n)。我们通过重复平方结果来实现这一点，因为# (result**2)**(2**(n-1)) = result**(2**n)而n > 1: result = result*result n- = 1 #最后，我们处理负x的情况。if(x < 0)：返回1/result else：返回结果
一个测试运行证实了这比exp2要好得多：
Math.exp (3.3)= 27.112638920657883 exp2(3.3)= 27.11262722548732 exp3(3.3)= 27.112638920657805
创建和使用模块
当程序变得更大、更复杂时，我们会想要开始将有用的函数组织到模块中，以便导入到程序中。在下一个例子中，我将编写一组处理素数的有用函数，然后将这些函数导入到一个简单的测试程序中。
首先，这是包含有用函数定义的模块：
defisprime (n,known_primes): "“”测试n是否为质数。第二个参数必须是一个质数列表，其中包括赢博体育小于根号n的质数。“”“对于known_primes中的p：如果p*p > n：返回True如果n % p == 0：返回False返回True def primeFactors(n,primes): ”“”计算并返回n的素数因子列表。"“” factors = [] x = n对于素数中的p：而x % p == 0: factors.append(p) x = x//p如果x >: 1: factors.append(x)返回factors deffindprimes （m,n，素数）：“”"构造并返回一个小于n的赢博体育素数的列表。素数必须是小于或等于m的素数的种子列表。myprime = primes[:]对于范围（m,n）中的num：如果isPrime（num, myprime）而不是myprime中的num: myprime .append（num）返回myprime
下面是使用该模块中函数的测试程序。
from prime_stuff import * p = [2,3,5,7] allPrimes = findPrimes(10,100,p) x = input（“输入2到10000之间的整数：”）如果isPrime(int(x),allPrimes): print（x +“是素数”）否则：factors = primeFactors(int(x),allPrimes) print（x +“不是素数”）print（“它的因子是”+ str（因子））
扩展的例子-数值积分
作为几何问题的积分
在微积分1的最后，我们介绍了函数f(x)的定积分。
这个量可以在几何上定义为面积：函数f(x)和x轴在x值a≤x≤b范围内的面积。正如你在微积分1中看到的，你可以使用几何技术来估计这个面积。最简单的几何技巧是将从a到b的区间划分为N + 1个等距点xi = a + h i，其中步长h等于(b-a)/N。
这个点的网格将区间[a，b]划分为子区间[xi， xi+1]，对于i = 0到n，我们可以在每个子区间上放置一个矩形，该矩形近似于子区间[xi， xi+1]上f(x)下面的面积。得到一个近似矩形的一种方法是给区间[xi， xi+1]上的矩形一个高度f（xi）。下图说明了这种设置。
把赢博体育这些近似矩形的面积加起来，我们就能得到我们想要的面积的粗略近似。
计算这个面积的程序
使用我上面描述的方法编写一个Python程序来估算面积相对容易。下面您将看到一个简单程序的代码，该程序使用这种方法来估计面积
这是一个很好的例子，因为在这个例子中曲线下面的面积是一个半径为2的四分之一圆，它的面积是π。因为我们可以预测确切的面积，我们可以在每次做矩形和时构造一个误差估计。程序重新计算不同N值的矩形面积之和，因此我们可以看到随着N的增加，误差缩小的速度有多快。
deff (x)：返回math.sqrt(4.0 - (2-x)*(2-x)) #通过使用步长为(b-a)/N的矩形#计算f(x)从a到b的积分#的近似值defrectanglearea (a,b，N): h = (b-a)/N；Sum = 0.0；对于范围(N)中的k: term = h*(f(a + k*h)) sum += term返回sum print（“N估计误差”）；N = 10 power = 1当power <= 6时：estimate = rectangleArea(0.0,2.0，N) error = math.fabs（math. fabs）pi -估计)print（"{:>10d} {:f} {:.16f}".format（N，估计，错误））N = N*10次方+= 1
下面是这个程序产生的结果。
N估计误差10 2.904518 0.2370743273414746 100 3.120417 0.0211756218107477 1000 3.139555 0.0020371866787654 10000 3.141391 0.0002011759784653 100000 3.141573 0.0000200371878285 1000000 3.141591 0.0000020011759911
正如您可以清楚地看到的，我们每增加10个细分，错误的大小就会减少10倍。这些结果的一个有点令人失望的方面是，误差降低得不够快。为了使误差减小到10-12，我们需要1012个子区间。计算赢博体育这些子区间的总面积所需的时间会非常大。
更复杂的方法
为了改善我们的误差估计，我们需要更好地将我们估计的面积与曲线下的实际面积相匹配。一种方法是用梯形代替矩形。在每个子区间[xi， xi+1]上，我们放置一个梯形，其左高由f（xi）给出，右高由f（xi+1）给出。下面的图片说明了梯形比矩形在匹配曲线下的面积方面做得更好。
一个这样的梯形的面积是h (f(xi) + f(xi+1))/2。
重写上面所示的面积计算程序以使用此公式进行面积估计是一件相对简单的事情。这是产生的结果。
N估计误差10 3.104518 0.0370743273414744 100 3.140417 0.0011756218107473 1000 3.141555 0.0000371866787638 10000 3.141591 0.0000011759784582 100000 3.141593 0.0000000371877436 1000000 3.141593 0.00000000371877436 0.0000003.141593 0.0000000011760504
这好多了，但还不够完美。Python浮点数可以表示精度为16位的数字，但我们只能将结果精确到10位左右。把结果的最后几位数挤出来就更难了。一种方法是继续增加细分的数量。这是不实际的，因为计算表中的最后一行需要一秒钟以上的时间，而且看起来我们需要使用100倍的细分才能得到精确到16位的结果。
理查森外推法
上面的情况说明了计算科学中的一个经典情况。一个相对粗糙的方法往往可以通过一些简单的调整得到显著的改进。然而，要达到较高的性能标准，就必须部署更多特殊的、专门的技术。这些笔记的其余部分将带你通过这种技术的发展，称为理查德森外推法。通过使用这种技术，我们将能够以更少的计算获得更高的性能。
理查德森外推从考虑梯形规则中的误差项开始。当我们使用梯形规则来估计N个子区间和步长为h的积分时，梯形规则通常会产生精确到与h2成正比的因子的结果。
这里用的符号O（h2）表示误差是h2阶的。另一种写法是把误差项写成h的不同次幂的组合
（这有点超出了这里讨论的范围，但事实证明，梯形法则的误差项只包含h的偶次幂。）
下一步是运用理查森外推法，这是一种奇怪但有效的代数技巧。在这个技巧中，我们把误差关系写了两次，一次是针对N个梯形，就像我们上面看到的，第二次是针对2n个梯形。将梯形的数量增加2倍，将步长从h减小到h/2，因此我们得到误差的两个表达式：
第二个方程乘以4，然后减去第一个方程，得到
我们可以把它写成
或
通过使用这个简单的代数技巧，我们已经从一个误差项大小为0 （h2）的方法切换到一个误差项大小为0 （h4）的方法。
更好的是，这个技巧可以反复使用。我就不赘述细节了，但是迭代这个技巧是一种叫做Romberg积分的积分技术的基础。下面是Romberg技术的总结。函数R（k,j）表示该技术的第j次迭代，从N = 2k-1步的梯形规则开始。
R(a,b,k,1) =梯形面积（a,b,2k-1）
下面是一个Python程序，它使用这个递归公式来计算积分估计。
deff (x)：返回math.sqrt(4.0 - (2-x)*(2-x)) #通过使用步长为(b-a)/N的梯形#计算f(x)从a到b的积分#的近似值defapezoidarea (a,b，N): h = (b-a)/N；Sum = 0.0；for k in range(N): term = h*(f(a + k*h)+f(a + (k+1)*h))/2 sum += term返回sum #递归定义的Romberg公式def Romberg (a,b,k,j)：如果j == 1：返回teapezoidarea (a,b,2**(k-1)) elif k < j：返回0否则：返回Romberg (a,b,k,j-1) + (Romberg (a,b,k,j-1))/(4**(j-1)-1) print（“ j估计误差”）；对于范围（2,6）中的j: estimate = romberg(0.0,2.0,4*j,j) error = math.fabs（math. fabs）。π-估计)打印(“{:> 3 d} {: f} {: .16f}”.format (j、估计、错误))
要解开这个程序正在做的事情需要一点努力。底部的循环是计算R（a,b, 4j,j）对于j在2到5的范围内。对于给定的j值，该方法从使用2k-1 = 24 j -1步的梯形规则估计开始。当j = 5时，这意味着使用220 - 1 = 524288步。这看起来很多，但请记住，在上一个示例程序中，我们将步数推到了一百万。
下面是这个程序返回的结果
J估计误差2 3.141275 0.0003172346542017 3 3.141588 0.0000043511837551 4 3.141593 0.0000000660130737 5 3.141593 0.0000000010239907
这个程序还有最后一个缺陷，这个缺陷严重地阻碍了从j= 4开始的计算。这里的问题是过度递归。从上面的代码中可以看到，对romberg函数的一个典型调用会触发三个额外的递归调用。如果这些依次触发对romberg函数的更多调用，那么最终可能会对该函数进行大量调用。
为了解决这个问题，我们可以使用一个简单的缓存策略。每次我们为k和j的特定组合计算romberg（a,b,k,j）的值时，我们应该记录下这个结果。下次我们需要这个值时，我们可以简单地查找结果，而不必重新计算任何东西。
程序的最终版本使用了这种缓存策略。我们建立了一个二维列表R，它可以存储k和j的多种组合的值。romberg函数然后将大部分工作交给递归内部函数r_int。r_int通过首先检查缓存来进行优化计算，看看k和j的组合是否已经记录了结果。如果是，它立即停止并返回缓存的值。如果没有，则返回递归公式来计算结果，然后在返回答案之前缓存该结果。
deff (x)：返回math.sqrt(4.0 - (2-x)*(2-x)) #通过使用步长为(b-a)/N的梯形#计算f(x)从a到b的积分#的近似值defapezoidarea (a,b，N): h = (b-a)/N；Sum = 0.0；对于k在range(N)内：term = h*(f(a + k*h)+f(a + (k+1)*h))/2 sum += term return sum #递归定义的Romberg公式与缓存# r_int将计算R[k，j]，但记住# result在R列表中以加快后续计算。def r_int(a,b,k,j，R): if R[k][j] != -1: #我们以前见过这个k，j吗？返回R [k] [j] #如果是,就返回缓存的值#如果没有,计算,缓存,并返回结果如果j = = 1: R [k] [j] = trapezoidArea (a, b, 2 * * (k - 1)) elif k <珍:R [k] [j] = 0: R [k] [j] = r_int (a, b, k, j - 1, R) + (r_int (a, b, k, j - 1, R) -r_int (a, b, k - 1、j - 1 R)) / (4 * * (j - 1) 1)返回R [k] [j] . def伯格(a, b, k, j): R = [[1 c的范围(j + 1)] R的范围(k + 1)]返回r_int (a, b, k, j, R)打印(“j估计错误”);对于range（3,22,2）中的j: estimate = romberg(0.0,2.0,j,j) error = math.fabs（math. fabs）。π-估计)打印(“{:> 3 d} {: f} {: .16f}”.format (j、估计、错误))
注意，由于我们现在有了一个更有效的方法，我们可以安全地扩展我们计算的j的范围。下面是这个程序返回的结果。
J估计误差3 3.090764 0.0508290045413791 5 3.135506 0.0060864702273435 7 3.140835 0.0007579750725255 9 3.141498 0.0000946982204093 11 3.141581 0.0000118360657910 13 3.141591 0.0000014794727754 15 3.141592 0.0000001849330067 17 3.141593 0.0000000231165633 19 3.141593 0.0000000028895282 21 3.141593 0.0000000003611502