第3 章 一些简单的数值程序

我们已经学习了一些Python语言基础结构,现在使用这些结构编写一些简单的程序。通过这
种方式,我们再顺便介绍几种语言结构和算法技术。
3.1 穷举法
当一个整数存在整数立方根时,图3-1中的代码会对其进行输出。如果输入不是一个完全立
方数,则输出一个消息进行说明。

#寻找完全立方数的立方根
x = int(input('Enter an integer: '))
ans = 0
while ans**3 < abs(x):
ans = ans + 1
if ans**3 != abs(x):
print(x, 'is not a perfect cube')
else:
if x < 0:
ans = -ans
print('Cube root of', x,'is', ans)

图3-1 使用穷举法求立方根

那么,对于何种x值,程序能正常结束呢?答案是“所有整数”。证明方法非常简单。
 表达式ans**3的值从0开始,并随着每次循环逐渐变大;
 当这个值达到或超过abs(x)时,循环结束;
 因为abs(x)的值总为正,所以循环结束前进行的迭代次数必然是有限的。
编写循环时,应该使用一个合适的递减函数。这个函数具有如下属性:
 它可以将一组程序变量映射为一个整数;
 进入循环时,它的值是非负的;
 当它的值≤0时,循环结束;
 每次循环它的值都会减小。

图3-1的while循环中,递减函数是什么呢?是abs(x) – ans**3。
下面,我们制造一些错误,看看会发生什么。首先,将语句ans = 0注释掉。Python解释器
会输出一条错误消息:
NameError: name 'ans' is not defined
因为解释器将ans绑定到任何对象之前,都要先找到与ans绑定的值。下面,我们还原ans的初始
化语句,将语句ans = ans + 1替换为ans = ans,再试着求8的立方根。如果你厌倦了漫长的等
待,可以按ctrl+c(同时按住ctrl键和c键),这样就可以回到shell的用户提示符。
下面,在循环开始部分添加一条语句:
print('Value of the decrementing function abs(x) - ans**3 is',
abs(x) - ans**3)
然后重新运行程序。这次程序会一次又一次地输出:
Value of the decrementing function abs(x) - ans**3 is 8
程序会永远运行下去,因为循环体没有减少ans**3和abs(x)之间的差距。遇到这种程序不
会正常结束的情况时,经验丰富的程序员经常会插入一些print语句,就像这次一样,测试递减
函数是否真的递减。
这个程序使用的算法技术称为穷举法,是猜测与检验算法的一个变种。我们枚举所有可能
性,直至得到正确答案或者尝试完所有值。乍看上去,这是一种极其愚蠢的解决方法。但令人
惊奇的是,穷举法经常是解决问题的最实用的方法。它实现起来特别容易,并且易于理解。还
有,在很多情况下,它的运行速度也足够快。请一定记得将你添加的pirnt语句删除或者注释
掉,并插入语句ans = ans + 1,然后试着求出1 957 816 251的立方根。程序几乎瞬间结束。
然后,再试试7 406 961 012 236 344 616。
眼见为实,即使需要几百万次猜测,也不是什么问题。现代计算机的速度真是太快了,它执
行一条指令只需1纳秒——10亿分之1秒。要想描述它有多快还真有些困难,光传输1英尺(约0.3
米)只需要1纳秒多一点。另外一种形容方式是,在你的声音传输100英尺的时间内,现代计算机
可以执行几百万条指令。
只是为了好玩,试着运行以下代码:

maxVal = int(input('Enter a postive integer: '))
i = 0
while i < maxVal:
i = i + 1
print(i)
看看你需要输入一个多大的整数,才能感受到在输出结果之前有个明显的时间间隔。
实际练习:编写一个程序,要求用户输入一个整数,然后输出两个整数root和pwr,满足0 <
pwr < 6,并且root**pwr等于用户输入的整数。如果不存在这样一对整数,则输出一条消息进
行说明。

3.2 for 循环
迄今为止,我们使用的while循环是高度程式化的,即都按照一个整数序列进行迭代。Python
提供了一种语言机制简化使用这种迭代方式的程序,这就是for循环。
for语句的一般形式如下(回忆一下,斜体文本是对程序中该处代码的一种描述,并不是实
际的代码):
for variable in sequence:
code block
for后面的变量被绑定到序列中的第一个值,并执行下面的代码块。然后变量被赋给序列中的第
二个值,再次执行代码块。该过程一直继续,直到穷尽这个序列或者执行到代码块中的break
语句。
绑定到变量的序列值通常使用内置函数range生成,它会返回一系列整数。range函数接受3
个整数参数:start、stop和step。生成一个数列:start、start + step、start + 2*step,
等等。如果step是正数,最后一个元素就是小于stop的最大整数start + i * step。如果step
是负数,最后一个元素就是大于stop的最小整数start + i * step。例如,表达式range(5, 40, 10)
会得到序列5, 15, 25, 35,range(40, 5, -10)会得到序列40, 30, 20, 10。如果省略第一个
参数,它会取默认值0;如果省略最后一个参数(步长) ,它会取默认值1。例如,range(0, 3)
和range(3)都会生成序列0, 1, 2。数列中的数值是以“按需产生”的原则生成的,所以即使
range(1000000)这样的表达式也只占用很少内存。①5.2节将更加深入地讨论range函数。
我们还可以通过字面量指定for循环中迭代的序列,如[0, 3, 2],但这种方式并不常用。
看下面的代码:
x = 4
for i in range(0, x):
print(i)
会输出:
0
1
2
3
再看看这段代码:
x = 4
for i in range(0, x):
print(i)
x = 5
它会引起一个问题,如果在循环中改变x的值,能否影响迭代次数?答案是“不能”。在for循环那行代码中,range函数的参数在循环的第一次迭代之前就已经被解释器求值,随后的迭代中不
会再次求值。
为了看一下实际运行情况,看下面的代码:
x = 4
for j in range(x):
for i in range(x):
print(i)
x = 2
会输出:
0
1
2
3
0
1
0
1
0
1
因为外层循环中的range函数只被求值一次,而内层循环中的range函数则在每次执行内层
for语句时都被求值。
图3-2中的代码重新实现了求立方根的穷举法。for循环中有一个break语句,使循环遍历完
成迭代序列中的所有元素之前终止。

#寻找完全立方数的立方根
x = int(input('Enter an integer: '))
for ans in range(0, abs(x)+1):
if ans**3 >= abs(x):
break
if ans**3 != abs(x):
print(x, 'is not a perfect cube')
else:
if x < 0:
ans = -ans
print('Cube root of', x,'is', ans)

图3-2 使用for循环和break语句

可以使用in操作符配合for循环语句,非常方便地遍历字符串中的字符。例如:
total = 0
for c in '12345678':
total = total + int(c)
print(total)
这段代码对字符串'12345678'中的数字进行求和,并输出最后的总数。

实际练习:假设s是包含多个小数的字符串,由逗号隔开,如s = '1.23, 2.4, 3.123'。编
写一个程序,输出s中所有数值的和。
3.3 近似解和二分查找
如果有人请你编写一个程序,求任意非负数的平方根,你应该怎么做?
你可能会要求一个更为精确的问题定义。例如,如果要找出2的平方根,程序应该怎么做?2
的平方根不是一个有理数,这意味着我们不可能将它的值表示成一个有限的数字序列(或一个
float类型的数),所以这个问题从一开始就是无解的。
实际上,我们想要的是一个能够找出近似解平方根的程序。也就是说,能够找到足够接近实
际平方根的近似解即可。我们会在后面仔细讨论这个问题。眼下,我们先认为“足够接近”的意
思就是,近似解位于实际解附近的一个常数范围内,这个常数我们称为epsilon。
图3-3中的代码实现了求近似平方根的算法。它使用了一个我们以前没有介绍过的操作符+=。
赋值语句ans += step在语义上等同于稍显冗长的代码ans = ans + step。操作符-=和*=也是
如此。

x = 25
epsilon = 0.01
step = epsilon**2
numGuesses = 0
ans = 0.0
while abs(ans**2 - x) >= epsilon and ans <= x:
ans += step
numGuesses += 1
print('numGuesses =', numGuesses)
if abs(ans**2 - x) >= epsilon:
print('Failed on square root of', x)
else:
print(ans, 'is close to square root of', x)

图3-3 使用穷举法求近似平方根

我们又一次使用了穷举法。请注意,这种求平方根的方法和你在中学学过的用铅笔求平方根
的方法完全不同。使用计算机解决问题的最好方法通常与手工解决问题的方法大相径庭。
上面的代码运行后会输出:
numGuesses = 49990
4.999000000001688 is close to square root of 25
这个程序没有发现25是个完全平方数,没有输出5,我们是不是应该大失所望呢?不,程序
做了它应该做的。尽管输出5是挺好的,但与输出一个足够接近5的数相比,并没有好到哪儿去。
如果我们设x = 0.25,你认为会发生什么情况?程序会找到一个足够接近0.5的数吗?不,
它会输出:

numGuesses = 2501
Failed on square root of 0.25
穷举法是一种查找技术,只在被查找集合中包含答案时才有效。这个例子中,我们对0和x
之间的值进行枚举。当x在0和1之间时,x的平方根不在这个区间内。改正的方法是,修改while
循环第一行中and的第二个操作数,得到:
while abs(ans**2 - x) >= epsilon and ans*ans <= x:
下面思考一下,程序会运行多长时间。迭代的次数依赖于答案与0的距离以及步长。大致说
来,程序会执行while循环至多x/step次。
我们在较大的数上试验一下这段代码,如x = 123 456。程序会运行一段时间,然后输出:
numGuesses = 3513631
Failed on square root of 123456
发生了什么?肯定存在一个浮点数,在0.01的范围内近似于123 456的平方根。为什么程序没
有找到它?问题在于我们的步长太大,程序跳过了所有合适的答案。试着将步长设为
epsilon**3,再运行一下程序。程序最终肯定会找到一个合适的答案,但你未必会有耐心等到
它运行结束。
程序大概会进行多少次猜测呢?步长为0.000001,123 456的平方根大约为351.36。这意味着
程序要进行351 000 000次左右的猜测,才能找到一个满意的答案。我们可以从一个接近答案的数
开始猜测,这样能快一些,但前提是我们知道答案。
是时候通过其他方法解决这个问题了。我们要改弦更张,选择一个更好的算法,而不是微调
现在的算法。但在这之前,我们先来看一个问题,这个问题乍看上去与求平方根完全没有关系。
思考这样一个问题:如何在英语字典中找出由给定字母序列开头的单词?穷举法在理论上可
以解决这个问题。你可以从第一个单词开始,检查每个单词,直到找到以这些字母开头的单词,
或者找遍所有单词。如果字典中有n个单词,那么平均n/2次检查就可以找到这个单词。如果这个
单词不在字典中,就需要n次检查。当然,那些使用纸质(不是在线版)字典查找单词的人绝对
不会使用这种方法。
幸运的是,出版字典的人会不辞劳苦地将单词按照字典顺序排列。这就使我们可以将字典翻
到单词可能存在的那一页(例如,以字母m开头的单词,可能在字典的中间页附近)。如果单词
开头的字母在字典顺序上位于这页中第一个单词的前面,我们就往前找;如果单词开头字母在这
页中最后一个单词的后面,我们就往后找。否则,我们就检查这些字母能否和本页中的某个单词
相匹配。
下面,我们将同样的理念应用于求x的平方根这个问题。假设我们知道x的平方根的一个非常
好的近似解位于0和max之间,就可以利用数值的全序性。也就是说,对于任意两个不同的数n1
和n2,都有n1 < n2,或者n1 > n2。所以,我们可以认为x的平方根位于下面直线上的某处:
0___________________________________________________________________max
并开始在这个区间内查找。因为我们不需要知道从哪里开始查找,所以可以从中间开始:
0_______________________________猜测________________________________max

如果这不是正确答案(多数时候不是),那么就看看它是太大还是太小。如果太大,我们就
可以知道答案肯定位于左侧;如果太小,我们就知道答案肯定位于右侧。然后可以在更小的区间
上重复这个过程。图3-4给出了这种算法的实现和测试。

x = 25
epsilon = 0.01
numGuesses = 0
low = 0.0
high = max(1.0, x)
ans = (high + low)/2.0
while abs(ans**2 - x) >= epsilon:
print('low =', low, 'high =', high, 'ans =', ans)
numGuesses += 1
if ans**2 < x:
low = ans
else:
high = ans
ans = (high + low)/2.0
print('numGuesses =', numGuesses)
print(ans, 'is close to square root of', x)

图3-4 使用二分查找求近似平方根
 

运行上面的代码,会输出:
low = 0.0 high = 25 ans = 12.5
low = 0.0 high = 12.5 ans = 6.25
low = 0.0 high = 6.25 ans = 3.125
low = 3.125 high = 6.25 ans = 4.6875
low = 4.6875 high = 6.25 ans = 5.46875
low = 4.6875 high = 5.46875 ans = 5.078125
low = 4.6875 high = 5.078125 ans = 4.8828125
low = 4.8828125 high = 5.078125 ans = 4.98046875
low = 4.98046875 high = 5.078125 ans = 5.029296875
low = 4.98046875 high = 5.029296875 ans = 5.0048828125
low = 4.98046875 high = 5.0048828125 ans = 4.99267578125
low = 4.99267578125 high = 5.0048828125 ans = 4.998779296875
low = 4.998779296875 high = 5.0048828125 ans = 5.0018310546875
numGuesses = 13
5.00030517578125 is close to square root of 25
请注意,这段程序找出的答案与前面算法并不相同。结果非常好,因为这个答案也完全符合问题
的要求。
更重要的是,我们发现每经过一次循环迭代,待查找空间都缩小了一半。因为这种算法每一
步都将查找空间分为两部分,所以称为二分查找。二分查找是对前面算法的一个重大改进,之前
的算法只能在每次迭代后将查找空间缩小一部分。
我们再试试x = 123 456,这次程序只进行30次猜测就找到一个可以接受的答案。
那x = 123 456 789呢?只需45次猜测。

我们应该使用这种算法计算平方根,这没什么好说的。此外,将算法中的2改成3,我们就可
以计算一个非负数的近似立方根。第4章会介绍一种语言机制,使我们可以将代码功能扩展为计
算任意次方的根。
实际练习:图3-4中,如果语句x = 25被替换为x = -25,代码会如何运行?
实际练习:应该如何修改图3-4中的代码,才能求出一个数的立方根?这个数既可以是正数,
也可以是负数。(提示:修改low保证答案位于待查找区域。)
3.4 关于浮点数
很多时候,float类型的数值是实数的一个非常好的近似。但“很多时候”并不代表所有情
况,这个功能失效时会引起不可思议的后果。例如,试着运行以下代码:
x = 0.0
for i in range(10):
x = x + 0.1
if x == 1.0:
print(x, '= 1.0')
else:
print(x, 'is not 1.0')s
与多数人一样,你可能会对输出的结果大吃一惊:
0.9999999999999999 is not 1.0
为什么会先执行else从句呢?
要想弄清楚发生这种情况的原因,我们应该知道浮点数在计算机中是如何表示的。为了搞清
这一点,我们需要了解二进制数。
第一次学习十进制数(也就是基数为10的数)时,我们就知道任何一个十进制数都可以用数
字序列0123456789中的数字表示。最右边的数字是100位,向左进一位是101位,以此类推。例如,
十进制数字序列302表示3 × 100 + 0 × 10 + 2 × 1。长度为n的序列可以表示多少个不同的数呢?长
度为1的序列可以表示10个数字(0~9)中的任何一个;长度为2的序列可以表示100个数(0~99)。
一般来说,长度为n的序列可以表示10n个不同的数。
二进制数(基数为2的数)的原理也是一样的。二进制数也可以表示成一个数字序列,其中
不是0就是1。这些数字经常称为位。最右边的数字是20位,向左进一位是21位,以此类推。比如,
二进制数字序列101表示1 × 4 + 0 × 2 + 1 × 1 = 5。那么长度为n的序列可以表示出多少个不同的数
呢?2n个。
实际练习:二进制数10011等于十进制中的哪个数?
可能是因为多数人都有10根手指,所以我们喜欢使用十进制表示数值。然而,所有现代计算
机系统都使用二进制表示数值。这并不是因为计算机生来有2根手指,而是因为容易制造硬件开
关,也就是仅有2种状态(开或闭)的设备。计算机使用二进制表示法,而人类使用十进制表示法,这就会导致认知上的不一致。
在几乎所有现代编程语言中,非整数数值都使用浮点数表示。现在,我们先假设计算机内部
使用的是十进制表示法,要将一个数表示成一个整数对:有效数字和指数。例如,1.949可以表
示为数对(1949, 3),它代表1949×103的积。
有效数字的数量决定了数值能被表示的精度。例如,如果只有两位有效数字,那么就无法准
确表示1.949。它会被转换成1.949的某个近似值,在这里是1.9。这种近似值称为舍入值。
现代计算机使用二进制表示法,而不是十进制表示法。我们使用二进制表示有效数字和指数,
而不是十进制,并且使用2作为指数的底数,而不是10。例如,0.625(5/8)会表示成数对(101, 11),
因为5/8是二进制的0.101,11是3的二进制表示,所以数对(101, 11)代表5 × 23 = 5/8 = 0.625。
那Python中写作0.1的十进制分数1/10呢?若使用4位有效数字,最好的表示方式是(0011,
101),等于3/32,也就是0.09375。如果有5位有效的二进制数字,可以将0.1表示成(11001, 1000),
等于25/256,也就是0.09765625。那么,需要多少位有效数字才能使用浮点数准确表示0.1呢?需
要无穷位!不存在两个整数sig和exp,使sig × 2exp = 0.1。所以无论Python(或任何一种语言)使
用多少位有效数字表示浮点数,都只能表示0.1的一个近似值。在多数Python版本中,使用53位精
度表示浮点数,所以为保存十进制0.1而使用的有效数字为:
11001100110011001100110011001100110011001100110011001
它等于十进制中的:
0.1000000000000000055511151231257827021181583404541015625
非常接近1/10,但并不是1/10。
回到本节开始的那段神秘代码,为什么
x = 0.0
for i in range(10):
x = x + 0.1
if x == 1.0:
print(x, '= 1.0')
else:
print(x, 'is not 1.0')
会输出:
0.9999999999999999 is not 1.0
我们现在知道,测试条件x == 1.0产生的结果是False,因为x绑定的值不是确切的1.0。如
果在else从句的末尾加上print x == 10.0 * 0.1这行代码,会输出什么呢?它会输出False,
因为在循环迭代中,Python至少有一次使用了所有有效数字并做了舍入。这可不是小学老师教给
我们的内容,但将0.1相加10次真的不等于10乘以0.1的值。①
顺便说一下,如果对浮点数进行舍入操作,可以使用round函数。表达式round(x, numDigits)会返回一个浮点数,等于将x保留小数点后numDigits位的舍入值。例如,print round(2**0.5, 3)
会输出1.414,作为2的平方根的近似值。
实数和浮点数之间的区别真的很重要吗?谢天谢地,大多数时候没有什么问题。几乎没有这
种情况:1.0可以接受但0.9999999999999999却不行。但是,需要注意对相等关系的检验。我们
已经看到,使用==比较两个浮点数会产生不可思议的结果。更合适的做法是,看看两个浮点数是
否足够接近,而不是这两个数是否相等。例如,编写代码时,abs(x - y) < 0.0001就比x == y
更好。
另一个需要注意的问题是累积的舍入误差。多数时候不会出现问题,因为计算机中保存的数
值有时候比预期值大一点,有时候又小一点。但在某些程序中,误差可能会沿着同一个方向累积。
3.5 牛顿-拉弗森法
最常用的近似算法通常被认为出自艾萨克·牛顿之手,称为“牛顿法”,但有时也称为“牛
顿拉弗森法”①。可以用它求出很多函数的实数根,但我们只用它求单变量多项式的实数根。要
想将这个方法扩展到多变量多项式,需要数学和算法两方面的知识。
单变量(按照惯例,我们用x表示变量)多项式或者是0,或者是一个有限数目的非零单项式
的和,如3x2 + 2x + 3。每一项(如3x2)都由一个常数(项的系数,这里是3)乘以变量(这里为x)
的非负整数次方(这里为2次方)组成。每项中变量的指数称为这一项的次数。多项式的次数就
是各项中的最大次数。比如,3(0次)、2.5x + 2(1次)和3x2(2次)。与之相反,2/x和x0.5都不是
多项式。
如果p是个多项式,r是个实数,我们就可以用p(r)表示当x = r时多项式的值。多项式p的根就
是方程p = 0的解,也就是实数r,使得p(r) = 0。例如,“求24的近似平方根”这个问题可以用公式
表示为:找到一个x,使得x2  24 ≈ 0。
牛顿证明了一个定理:如果存在一个值guess是多项式p的根的近似值,那么guess 
p(guess)/p'(guess)就是一个更好的近似值②,其中p'是p的一次导数。
对于任意的常数k和任意的系数c,多项式cx2 + k的一次导数是2cx。例如,x2 – k的一次导数
是2x。因此,如果当前的猜测是y,那么可以选择y – (y2  k)/2y作为下一个猜测。这种方法称为逐
次逼近。图3-5中的代码展示了如何使用这种思想快速找出近似平方根。
 

#利用牛顿拉弗森法寻找平方根
#寻找x,满足x**2-24在epsilon和0之间
epsilon = 0.01
k = 24.0
guess = k/2.0
while abs(guess*guess - k) >= epsilon:
guess = guess - (((guess**2) - k)/(2*guess))
print('Square root of', k, 'is about', guess)

图3-5 牛顿拉弗森法

实际练习:在牛顿拉弗森法的实现中添加一些代码,跟踪求平方根所用的迭代次数。在
这段代码的基础上编写一个程序,比较牛顿拉弗森法和二分查找法的效率(你会发现牛顿拉
弗森法效率更高)。


 

  • 23
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
很高兴回答您的问题!以下是第3MATLAB数值计算的一些习题解析: 1. 编写一个 MATLAB 程序,计算方程 x^3 - 2x - 5 = 0 的根。 解析:可以使用 fzero 函数来计算非线性方程的根。具体代码如下: ```matlab f = @(x) x^3 - 2*x - 5; x0 = 0; % 设定初始值 x = fzero(f, x0); % 计算方程的根 disp(['方程的根为:', num2str(x)]); ``` 2. 使用 MATLAB 计算定积分 $\int_0^1 x^2 dx$。 解析:可以使用 quad 函数来计算定积分。具体代码如下: ```matlab f = @(x) x.^2; % 定义被积函数 a = 0; % 积分下限 b = 1; % 积分上限 I = quad(f, a, b); % 计算定积分 disp(['定积分的值为:', num2str(I)]); ``` 3. 编写一个 MATLAB 程序,计算矩阵的行列式。 解析:可以使用 det 函数来计算矩阵的行列式。具体代码如下: ```matlab A = [1 2 3; 4 5 6; 7 8 9]; % 定义矩阵 detA = det(A); % 计算矩阵的行列式 disp(['矩阵的行列式为:', num2str(detA)]); ``` 4. 编写一个 MATLAB 程序,求解线性方程组 Ax=b。 解析:可以使用 linsolve 函数来求解线性方程组。具体代码如下: ```matlab A = [1 2 3; 4 5 6; 7 8 9]; % 定义系数矩阵 b = [1; 2; 3]; % 定义常数向量 x = linsolve(A, b); % 求解线性方程组 disp(['方程组的解为:', num2str(x')]); ``` 5. 编写一个 MATLAB 程序,计算向量的范数。 解析:可以使用 norm 函数来计算向量的范数。具体代码如下: ```matlab v = [1; 2; 3]; % 定义向量 p = 2; % 范数的阶数 n = norm(v, p); % 计算向量的范数 disp(['向量的范数为:', num2str(n)]); ``` 希望这些习题解析能对您有所帮助!

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

___Y1

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值