（1）理解科学计算实质并掌握Python语言的科学计算应用

（2）掌握常用科学计算库

（3）熟练运用numpy及scipy、matplotlib等计算库资源 

【实验准备】

Python核心科学计算库的导入、配置并熟悉相关对象

【实验内容】

1. 实验练习：Numpy 应用包

NumPy 中除了可以使用 numpy.transpose 函数来对换数组的维度，还可以使用 T 属性。。

例如有个 m 行 n 列的矩阵，使用 t() 函数就能转换为 n 行 m 列的矩阵。

实例

import numpy as np a = np.arange(12).reshape(3,4) print ('原数组：') print (a) print ('\n') print ('转置数组：') print (a.T)

输出结果如下：

原数组：

[[ 0 1 2 3]

[ 4 5 6 7]

[ 8 9 10 11]]

转置数组：

[[ 0 4 8]

[ 1 5 9]

[ 2 6 10]

[ 3 7 11]]

matlib.empty()

matlib.empty() 函数返回一个新的矩阵，语法格式为：

numpy.matlib.empty(shape, dtype, order)

参数说明：

1. shape: 定义新矩阵形状的整数或整数元组
2. Dtype: 可选，数据类型
3. order: C（行序优先） 或者 F（列序优先）

实例

import numpy.matlib import numpy as np print (np.matlib.empty((2,2))) # 填充为随机数据

输出结果为：

[[-1.49166815e-154 -1.49166815e-154]

[ 2.17371491e-313 2.52720790e-212]]

numpy.matlib.zeros()

numpy.matlib.zeros() 函数创建一个以 0 填充的矩阵。

实例

import numpy.matlib import numpy as np print (np.matlib.zeros((2,2)))

输出结果为：

[[0. 0.]

[0. 0.]]

numpy.matlib.ones()

numpy.matlib.ones()函数创建一个以 1 填充的矩阵。

实例

import numpy.matlib import numpy as np print (np.matlib.ones((2,2)))

输出结果为：

[[1. 1.]

[1. 1.]]

numpy.matlib.eye()

numpy.matlib.eye() 函数返回一个矩阵，对角线元素为 1，其他位置为零。

numpy.matlib.eye(n, M,k, dtype)

参数说明：

1. n: 返回矩阵的行数
2. M: 返回矩阵的列数，默认为 n
3. k: 对角线的索引
4. dtype: 数据类型

实例

import numpy.matlib import numpy as np print (np.matlib.eye(n = 3, M = 4, k = 0, dtype = float))

输出结果为：

[[1. 0. 0. 0.]

[0. 1. 0. 0.]

[0. 0. 1. 0.]]

numpy.matlib.identity()

numpy.matlib.identity() 函数返回给定大小的单位矩阵。

单位矩阵是个方阵，从左上角到右下角的对角线（称为主对角线）上的元素均为 1，除此以外全都为 0。

实例

import numpy.matlib import numpy as np # 大小为 5，类型位浮点型 print (np.matlib.identity(5, dtype = float))

输出结果为：

[[ 1. 0. 0. 0. 0.]

[ 0. 1. 0. 0. 0.]

[ 0. 0. 1. 0. 0.]

[ 0. 0. 0. 1. 0.]

[ 0. 0. 0. 0. 1.]]

numpy.matlib.rand()

numpy.matlib.rand() 函数创建一个给定大小的矩阵，数据是随机填充的。

实例

import numpy.matlib import numpy as np print (np.matlib.rand(3,3))

输出结果为：

[[0.23966718 0.16147628 0.14162 ]

[0.28379085 0.59934741 0.62985825]

[0.99527238 0.11137883 0.41105367]]

矩阵总是二维的，而 ndarray 是一个 n 维数组。 两个对象都是可互换的。

实例

import numpy.matlib import numpy as np i = np.matrix('1,2;3,4') print (i)

输出结果为：

[[1 2]

[3 4]]

实例

import numpy.matlib import numpy as np j = np.asarray(i) print (j)

输出结果为：

[[1 2]

[3 4]]

实例

import numpy.matlib import numpy as np k = np.asmatrix (j) print (k)

输出结果为：

[[1 2]

[3 4]]

2. 实验练习：Numpy 切片和索引

ndarray对象的内容可以通过索引或切片来访问和修改，与 Python 中 list 的切片操作一样。

ndarray 数组可以基于 0 - n 的下标进行索引，切片对象可以通过内置的 slice 函数，并设置 start, stop 及 step 参数进行，从原数组中切割出一个新数组。

实例

import numpy as np a = np.arange(10) s = slice(2,7,2) # 从索引 2 开始到索引 7 停止，间隔为2 print (a[s])

输出结果为：

[2 4 6]

以上实例中，我们首先通过 arange() 函数创建 ndarray 对象。 然后，分别设置起始，终止和步长的参数为 2，7 和 2。

我们也可以通过冒号分隔切片参数 start:stop:step 来进行切片操作：

实例

import numpy as np a = np.arange(10) b = a[2:7:2] # 从索引 2 开始到索引 7 停止，间隔为 2 print(b)

输出结果为：

[2 4 6]

冒号 : 的解释：如果只放置一个参数，如 [2]，将返回与该索引相对应的单个元素。如果为 [2:]，表示从该索引开始以后的所有项都将被提取。如果使用了两个参数，如 [2:7]，那么则提取两个索引(不包括停止索引)之间的项。

实例

import numpy as np a = np.arange(10) # [0 1 2 3 4 5 6 7 8 9] b = a[5] print(b)

输出结果为：

5

实例

import numpy as np a = np.arange(10) print(a[2:])

输出结果为：

[2 3 4 5 6 7 8 9]

实例

import numpy as np a = np.arange(10) # [0 1 2 3 4 5 6 7 8 9] print(a[2:5])

输出结果为：

[2 3 4]

多维数组同样适用上述索引提取方法：

实例

import numpy as np a = np.array([[1,2,3],[3,4,5],[4,5,6]]) print(a) # 从某个索引处开始切割 print('从数组索引 a[1:] 处开始切割') print(a[1:])

输出结果为：

[[1 2 3]

[3 4 5]

[4 5 6]]

从数组索引 a[1:] 处开始切割

[[3 4 5]

[4 5 6]]

切片还可以包括省略号 …，来使选择元组的长度与数组的维度相同。 如果在行位置使用省略号，它将返回包含行中元素的 ndarray。

实例

import numpy as np a = np.array([[1,2,3],[3,4,5],[4,5,6]]) print (a[...,1]) # 第2列元素 print (a[1,...]) # 第2行元素 print (a[...,1:]) # 第2列及剩下的所有元素

输出结果为：

[2 4 5]

[3 4 5]

[[2 3]

[4 5]

[5 6]]

3. scipy库应用：scipy.cluster 聚类scipy.constants 数学常量scipy.fft 快速傅里叶变换scipy.integrate 积分scipy.interpolate 插值scipy.io 数据输入输出scipy.linalg 线性代数scipy.misc 图像处理scipy.ndimage N 维图像scipy.odr 正交距离回归scipy.optimize 优化算法scipy.signal 信号处理scipy.sparse 稀疏矩阵scipy.spatial 空间数据结构和算法scipy.special 特殊数学函数scipy/stats 统计函数

（以上为scipy子库，请同学熟悉库名）

实例：求解f(x)=2*sin(x)-x+1参考程序：

def f1(x):

       return np.sin(x)*2-x+1

value = fsolve(f1,[2])

print(value)

4.Numpy实例 卫星GPS定位

 上述6个卫星S1-S6和地球在高速运动，从卫星发出的位置信息以光速传输到GPS接收端需要一定的时间。  假设(x,y,z,t)表示R当前的位置, t是R的相对时间，卫星S1（发出信号时刻）到（当前接收时刻)满足以下关系(其中c是光速)。 (x-3)^2 + (y-2)^2 + (z-3)^2 = [(10010.00692286 – t)*c]^2， 该公式表示以(x, y, z,t)为参数的(欧式空间距离)与信号传输距离相等。对于卫星S1,S2,…,S6，满足方程组:...(1)

其中，光速为常数c=0.299792458km/us，上述方程组是非线性的，但很容易将所有二次项都消去（每个公式减去第一个公式），从而得到：

此时，上述等式变成了A*X=B形式，根据线性代数方法，X=A-1*B，即只需对系数矩阵求逆，再乘以常数矩阵便可以得到方程组的解。GPS定位的问题建模： 上面给出了GPS的定位原理，如何利用计算机辅助GPS的定位计算呢？ 以6颗卫星为例，GPS定位计算问题的IPO模式----描述如下： 输入：6颗卫星的欧式坐标和信号时间戳 处理：GPS定位算法 输出：GPS接收设备的地理坐标和当前时间假设第i颗卫星的坐标和时间戳表示为(x_i, y_i ,z_i ,t_i )，结合上述例子，GPS定位算法可以描述为如下公式：

下面将使用Numpy函数库实现上述矩阵操作。程序中用到的函数：

numpy.dot(a,b)：计算矩阵a与矩阵b的点积

numpy.linalg.inv(a)：求矩阵a 的逆矩阵

GPS定位的程序实现 Python代码如下： 其中zeros是NumPy提供的函数，用来建立指定维度的数组, zeros用来生成数组x用来存储接受来自外部输入的六颗卫星坐标， 数组a，b用来存放前面算法中的系数矩阵，

参考程序：【请阅读并分析写出主要代码】

from numpy import *

def main_GPSLocation():

    i = 1

    c = 0.299792458  # 光速 0.299792458km/us

    x = zeros((6, 4)) #存储6个卫星的（x,y,z,t）参数

    while i<=6:

        print(" %s %d" % ("please input (x,y,z,t) of group",i) )

        temp=input()

        x[i-1]=temp.split()

        j=0

        while j<4:

            x[i-1][j]=float(x[i-1][j])

            j=j+1

        i=i+1

    a=zeros((4,4)) #系数矩阵

    b=zeros((4,1)) #常数项

    j=0

    while j<4:

        a[j][0]=2*(x[5][0]-x[j][0])

        a[j][1]=2*(x[5][1]-x[j][1])

        a[j][2]=2*(x[5][2]-x[j][2])

        a[j][3]=2*c*c*(x[j][3]-x[5][3])

        b[j][0]=x[5][0] * x[5][0] - x[j][0] * x[j][0] + \

                x[5][1] * x[5][1] - x[j][1] * x[j][1] + \

                x[5][2] * x[5][2] - x[j][2] * x[j][2] + \

            c*c*(x[j][3] * x[j][3] - x[5][3] * x[5][3])

        j=j+1

    a_ni=linalg.inv(a) #系数矩阵求逆

    print(dot(a_ni,b))

main_GPSLocation()

【思考题】

1. 请思考Python numpy库中针对数组的操作与列表、元组等序列数据处理有何联系及区别？
2. 如果想在scipy子库linalg中增加一些针对矩阵的处理代码，如何实现？请谈谈自己的。

Python的NumPy库提供了对数组进行高效处理的功能，与列表和元组等序列数据处理有着联系和区别。

联系：

• 数据处理：NumPy的数组（numpy.array）和Python的列表（list）都可以存储多个元素，并进行索引、切片等操作。
• 元素访问：两者都可以通过索引访问元素，但NumPy数组支持更多的高效操作，例如支持向量化操作和广播。
• 数据类型：NumPy数组可以指定特定的数据类型（例如int32、float64等），而列表则可以包含不同类型的对象。

区别：

• 性能：NumPy数组在执行数值计算时通常比Python列表更高效，因为NumPy是基于C语言编写的，拥有更多的优化和内置函数。
• 维度：NumPy数组可以是多维的（例如二维矩阵），而Python列表通常是一维的，可以嵌套以实现多维效果。
• 功能：NumPy提供了许多数学和统计函数，例如求和、平均值、矩阵乘法等，这些函数对于数组操作非常方便。

关于在SciPy子库（比如scipy.linalg）中增加针对矩阵处理的代码，一般可以通过以下步骤实现：

1. 了解子库的结构和需求：首先需要深入了解scipy.linalg子库中已有的功能和提供的矩阵处理方法，确保你的添加是有意义且没有重复的。

2. 贡献代码：你可以参考SciPy的贡献指南，一般来说，你可以先在GitHub上找到SciPy仓库，进行Fork，然后修改代码并提交Pull Request。在修改代码之前，最好先在邮件列表或者GitHub上提出你的想法，以便其他开发者提供反馈和指导。

3. 编写文档和测试：编写清晰的文档说明你添加的新功能，包括用法、示例和可能的注意事项。同时，编写测试用例来验证新功能的正确性和稳定性。

4. 遵循社区规范：在贡献代码时，遵循SciPy项目的代码风格、规范和最佳实践，确保你的代码与现有代码风格一致，并通过项目的审核和测试流程。