Growing

日拱一卒。

发表于 更新于
本文字数： 111k 阅读时长 ≈ 1:41

此篇记录在python实操过程中遇见的问题以及最终的解决办法

1
2
3
4
5
6
7
8
9
# 导入常用的包

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import scipy as si

%matplotlib inline

写在最前

  • 善用帮助：函数和文档
  • 多用函数和循环，减少代码量

读写文件的几种方式

  1. python

    • f.read();f.readline();f.readlines()
    • f.write();f.writelines();f.close()
    • f.seek();f.tell();

  2. numpy

    • numpy.loadtxt();numpy.savetxt()
    • numpy.tofile();numpy.fromfile()
    • numpy.genfromtxt();numpy.recfromcsv()
    • numpy.save();numpy.load() .npy
    • numpy.savez();numpy.loadz() .npz
    • numpy.savez_compressed()

  3. scipy

    • scipy.io.savemat();scipy.io.loadmat()

  4. pandas

    • pd.read_csv()
    • pd.read_hdf()
    • pd.read_excel()

读取图像

  1. scipy

    • scipy.misc.imread()

  2. matplotlib.pyplot

    • plt.imread()

  3. scikit-image:

    • skimage.io
1
# help(np.loadtxt)

速查表

  • Bokeh:Bokeh 是 Python 的交互式可视图库，用于生成在浏览器里显 示的大规模数据集高性能可视图。
  • jupyter notebook:Jupyter将代码与文本封装为三种类型的单元格：Markdown、代码与 NBConvert
  • Keras:Keras是强大、易用的深度学习库，基于Theano和TensorFlow提供 了高阶神经网络API，用于开发和评估深度学习模型。
  • Matplotlib:Matplotlib 是 Python 的二维绘图库，用于生成符合出版质量或 跨平台交互环境的各类图形
  • Numpy:Numpy 是 Python 数据科学计算的核心库，提供了高性能的多维数组对象及处 理数组的工具。
  • Pandas:Pandas 是基于 Numpy 创建的 Python 库，为 Python 提供 了易于使用的数据结构和数据分析工具。
  • Skit-Learn:Scikit-learn 是开源的 Python 库，通过统一的界面实现 机器学习、预处理、交叉验证及可视化算法。
  • Seaborn:Seaborn 是基于 matplotlib 开发的高阶Python 数据可视图库， 用于绘制优雅、美观的统计图形。
  • 导入数据：大多数情况下，都是用 Numpy 或 Pandas 导入数据。

numpy

numpy 索引

  • 切片和索引
  • 高级索引
    • 整数数组索引
    • 布尔索引
    • 花式索引
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
import numpy as np
 
a = np.arange(10)
s = slice(2,7,2)   # 从索引 2 开始到索引 7 停止，间隔为2
print (a[s])
[2 4 6]

 
a = np.arange(10)
print(a[2:])
[2 3 4 5 6 7 8 9]


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


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]]


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

[[1 2 3]
 [3 4 5]
 [4 5 6]] 

[2 4 5] 

[3 4 5] 

[[2 3]
 [4 5]
 [5 6]]

 
x = np.array([[1,  2],  [3,  4],  [5,  6]])
print(x,'\n')
y = x[[0,1,2],  [0,1,0]]  
print (y)

[[1 2]
 [3 4]
 [5 6]] 

[1 4 5]

 
x = np.array([[  0,  1,  2],[  3,  4,  5],[  6,  7,  8],[  9,  10,  11]])  
print ('我们的数组是：' )
print (x)
print ('\n')
rows = np.array([[0,0],[3,3]]) 
cols = np.array([[0,2],[0,2]]) 
y = x[rows,cols]  
print  ('这个数组的四个角元素是：')
print (y)

我们的数组是：
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]


这个数组的四个角元素是：
[[ 0  2]
 [ 9 11]]

 
a = np.array([[1,2,3], [4,5,6],[7,8,9]])
b = a[1:3, 1:3]
c = a[1:3,[1,2]]
d = a[...,1:]
print(b)
print(c)
print(d)

[[5 6]
 [8 9]]
[[5 6]
 [8 9]]
[[2 3]
 [5 6]
 [8 9]]

 
x = np.array([[  0,  1,  2],[  3,  4,  5],[  6,  7,  8],[  9,  10,  11]])  
print ('我们的数组是：')
print (x)
print ('\n')
# 现在我们会打印出大于 5 的元素  
print  ('大于 5 的元素是：')
print (x[x >  5])

我们的数组是：
[[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]]


大于 5 的元素是：
[ 6  7  8  9 10 11]

 
a = np.array([np.nan,  1,2,np.nan,3,4,5]) 
print(a,'\n')
print (a[~np.isnan(a)])

[nan  1.  2. nan  3.  4.  5.] 

[1. 2. 3. 4. 5.]

 
a = np.array([1,  2+6j,  5,  3.5+5j]) 

print (a[np.iscomplex(a)])

[2. +6.j 3.5+5.j]

a = np.arange(3*4*4).reshape(3,4,4)
a

array([[[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11],
        [12, 13, 14, 15]],

       [[16, 17, 18, 19],
        [20, 21, 22, 23],
        [24, 25, 26, 27],
        [28, 29, 30, 31]],

       [[32, 33, 34, 35],
        [36, 37, 38, 39],
        [40, 41, 42, 43],
        [44, 45, 46, 47]]])
1
2
b = a[1]
b
1
2
3
4
array([[16, 17, 18, 19],
       [20, 21, 22, 23],
       [24, 25, 26, 27],
       [28, 29, 30, 31]])
1
2
b[1:3,1:3]=True
b
1
2
3
4
array([[16, 17, 18, 19],
       [20,  1,  1, 23],
       [24,  1,  1, 27],
       [28, 29, 30, 31]])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
mask = np.empty((256,256))


for i in range(1,257):
    for j in range(1,257):
        if ((i-127.5)**2 + (j-98)**2) <= (53.8**2):
            mask[i-1,j-1]=True
          
        else:
            mask[i-1,j-1]=False


import matplotlib.pyplot as plt


plt.imshow(mask)
1
<matplotlib.image.AxesImage at 0x17ab71641f0>

png

1
mask
1
2
3
4
5
6
7
array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])
1
2
temp = np.full((256,256),5)
plt.imshow(temp)
1
<matplotlib.image.AxesImage at 0x17ab718ddf0>

png

1
mask = mask.astype('bool')
1
mask
1
2
3
4
5
6
7
array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])
1
2
print(np.shape(temp[mask]))
print(np.pi*53.8**2)
1
2
(9092,)
9093.151440256439
1
plt.imshow(mask)
1
<matplotlib.image.AxesImage at 0x17755787100>

png

1

1

创建数组array和矩阵matrix

numpy中二者几乎等效，且array更高效

1
2
3
4
import numpy as np
a = np.array([2,3,4])
print(a)
print(a.shape,a.ndim,a.dtype,a.dtype.name,a.itemsize,a.size,a.data)
1
2
[2 3 4]
(3,) 1 int32 int32 4 3 <memory at 0x000002B1B0350C40>
1
2
3
b = np.arange(12).reshape(4,3)
b
# array, zeros, zeros_like, ones, ones_like, empty, empty_like, arange, linspace,
1
2
3
4
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])

基本运算

  • 广播运算
  • 元素乘法
  • 矩阵乘法
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 广播运算

# 对应元素进行运算：
a = np.array( [20,30,40,50] )
b = np.arange( 4 )
print(b)
c = a-b
print(c)
print(b**2)
print(b+5)

print(10*np.sin(a))

a<35
1
2
3
4
5
6
7
8
9
10
11
[0 1 2 3]
[20 29 38 47]
[0 1 4 9]
[5 6 7 8]
[ 9.12945251 -9.88031624  7.4511316  -2.62374854]





array([ True,  True, False, False])
1
2
# 对应元素相乘
a*b
1
array([  0,  30,  80, 150])
1
2
3
4
5
6
7
8
# 矩阵乘法
A = np.array( [[1,1],
               [0,1]] )
B = np.array( [[2,0],
               [3,4]] )
print(A * B)                      # elementwise product: 点乘
print(A @ B)               # matrix product
print(A.dot(B))                  # another matrix product
1
2
3
4
5
6
[[2 0]
 [0 4]]
[[5 4]
 [3 4]]
[[5 4]
 [3 4]]
1
2
temp = np.arange(4).reshape(2,2)
temp
1
2
array([[0, 1],
       [2, 3]])
1
temp.mean()
1
1.5
1
temp.mean(axis=0)
1
array([1., 2.])
1
2
temp.mean(axis=1)
# a.max,a.min,a.var,a.std,...
1
array([0.5, 2.5])
1
2
temp1 = np.arange(8).reshape(2,2,2)
temp1
1
2
3
4
5
array([[[0, 1],
        [2, 3]],

       [[4, 5],
        [6, 7]]])
1
temp1[0,0,0]
1
0
1
temp1[1,0,0]
1
4
1
2
temp1.mean(axis=0)
# 即沿着轴向取平均
1
2
array([[2., 3.],
       [4., 5.]])
1
temp1.mean(axis=1)
1
2
array([[1., 2.],
       [5., 6.]])
1
temp1.mean(axis=2)
1
2
array([[0.5, 2.5],
       [4.5, 6.5]])
1

1
# 广播运算
1

1

1

线性代数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import numpy as np
a = np.array([[1.0, 2.0], [3.0, 4.0]])
print('a=\n',a)

print("a.T=\n",a.transpose()) # 转置
print("a'=\n",np.linalg.inv(a)) # 求逆矩阵

u = np.eye(2) # unit 2x2 matrix; "eye" represents "I"
print('I=\n',u) # 单位矩阵
j = np.array([[0.0, -1.0], [1.0, 0.0]])
print('j=\n',j)
print(j @ j)   # matrix product

print(np.trace(u))  # trace :矩阵的迹
y = np.array([[5.], [7.]])
print('result=\n',np.linalg.solve(a, y)) # 求线性方程

print('eig=\n',np.linalg.eig(j)) # 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
a=
 [[1. 2.]
 [3. 4.]]
a.T=
 [[1. 3.]
 [2. 4.]]
a'=
 [[-2.   1. ]
 [ 1.5 -0.5]]
I=
 [[1. 0.]
 [0. 1.]]
j=
 [[ 0. -1.]
 [ 1.  0.]]
[[-1.  0.]
 [ 0. -1.]]
2.0
result=
 [[-3.]
 [ 4.]]
eig=
 (array([0.+1.j, 0.-1.j]), array([[0.70710678+0.j        , 0.70710678-0.j        ],
       [0.        -0.70710678j, 0.        +0.70710678j]]))

python图像处理

1.NumPy +SciPy
2.scikit-image

scipy 图像处理

1
2
3
4
5
6
from scipy import misc
import matplotlib.pyplot as plt

face = misc.face()
plt.imshow(face)
plt.show()

png

几何变换

1
2
3
4
5
6
7
8
9
10
11
from scipy import ndimage

face = misc.face(gray=True)
plt.imshow(face,cmap='gray')

shifted_face = ndimage.shift(face, (50, 50))
shifted_face2 = ndimage.shift(face, (50, 50), mode='nearest')
rotated_face = ndimage.rotate(face, 30)
cropped_face = face[50:-50, 50:-50]
zoomed_face = ndimage.zoom(face, 2)
zoomed_face.shape
1
(1536, 2048)

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
plt.figure(figsize=(10,12))
plt.subplot(321)  
plt.imshow(face, cmap='gray')  
plt.axis('off')

plt.subplot(322)  
plt.imshow(shifted_face, cmap=plt.cm.gray)  
plt.axis('off')

plt.subplot(323)  
plt.imshow(shifted_face2, cmap=plt.cm.gray)  
plt.axis('off')

plt.subplot(324)  
plt.imshow(rotated_face, cmap=plt.cm.gray)  
plt.axis('off')

plt.subplot(325)  
plt.imshow(cropped_face, cmap=plt.cm.gray)  
plt.axis('off')

plt.subplot(326)  
plt.imshow(zoomed_face, cmap=plt.cm.gray)  
plt.axis('off')
1
(-0.5, 2047.5, 1535.5, -0.5)

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
# Load some data
from scipy import misc
face = misc.face(gray=True)

# Apply a variety of transformations
from scipy import ndimage
from matplotlib import pyplot as plt
shifted_face = ndimage.shift(face, (50, 50))
shifted_face2 = ndimage.shift(face, (50, 50), mode='nearest')
rotated_face = ndimage.rotate(face, 30)
cropped_face = face[50:-50, 50:-50]
zoomed_face = ndimage.zoom(face, 2)
zoomed_face.shape

plt.figure(figsize=(15, 3))
plt.subplot(151)
plt.imshow(shifted_face, cmap=plt.cm.gray)
plt.axis('off')

plt.subplot(152)
plt.imshow(shifted_face2, cmap=plt.cm.gray)
plt.axis('off')

plt.subplot(153)
plt.imshow(rotated_face, cmap=plt.cm.gray)
plt.axis('off')

plt.subplot(154)
plt.imshow(cropped_face, cmap=plt.cm.gray)
plt.axis('off')

plt.subplot(155)
plt.imshow(zoomed_face, cmap=plt.cm.gray)
plt.axis('off')

plt.subplots_adjust(wspace=.05, left=.01, bottom=.01, right=.99, top=.99)

plt.show()

png

图像过滤

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
# Load some data
from scipy import misc
face = misc.face(gray=True)
face = face[:512, -512:]  # crop out square on right

# Apply a variety of filters
from scipy import ndimage
from scipy import signal
from matplotlib import pyplot as plt

import numpy as np
noisy_face = np.copy(face).astype(np.float)
noisy_face += face.std() * 0.5 * np.random.standard_normal(face.shape)
blurred_face = ndimage.gaussian_filter(noisy_face, sigma=3)
median_face = ndimage.median_filter(noisy_face, size=5)
wiener_face = signal.wiener(noisy_face, (5, 5))

plt.figure(figsize=(12, 3.5))
plt.subplot(141)
plt.imshow(noisy_face, cmap=plt.cm.gray)
plt.axis('off')
plt.title('noisy')

plt.subplot(142)
plt.imshow(blurred_face, cmap=plt.cm.gray)
plt.axis('off')
plt.title('Gaussian filter')

plt.subplot(143)
plt.imshow(median_face, cmap=plt.cm.gray)
plt.axis('off')
plt.title('median filter')

plt.subplot(144)
plt.imshow(wiener_face, cmap=plt.cm.gray)
plt.title('Wiener filter')
plt.axis('off')

plt.subplots_adjust(wspace=.05, left=.01, bottom=.01, right=.99, top=.99)

plt.show()

png

数学形态学运算

  1. 二值形态学
    1. Erosion 腐蚀； scipy.ndimage.binary_erosion()
    2. Dilation 膨胀； scipy.ndimage.binary_dilation()
    3. Opening 开运算； scipy.ndimage.binary_opening()
    4. Closing 闭运算； scipy.ndimage.binary_closing()
  2. 灰度形态学
1

图像处理：scikit-image

1
2
3
4
5
6
import numpy as np
check = np.zeros((9, 9))
check[::2, 1::2] = 1
check[1::2, ::2] = 1
import matplotlib.pyplot as plt
plt.imshow(check, cmap='gray', interpolation='nearest') 
1
<matplotlib.image.AxesImage at 0x17ab70e3af0>

png

1
2
import skimage
from skimage import data  # most functions are in subpackages
1
2
3
4
5
6
7
8
camera = data.camera()
camera.dtype
plt.imshow(camera)
camera.shape

from skimage import filters
filtered_camera = filters.gaussian(camera, 1)
type(filtered_camera)   
1
numpy.ndarray

png

图像处理：scipy.ndimage

1
2
3
4
5
6
7
from scipy import misc
f = misc.face()
misc.imsave('face.png', f) # uses the Image module (PIL)

import matplotlib.pyplot as plt
plt.imshow(f)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
---------------------------------------------------------------------------

AttributeError                            Traceback (most recent call last)

<ipython-input-8-d3ddc3b982b6> in <module>
      1 from scipy import misc
      2 f = misc.face()
----> 3 misc.imsave('face.png', f) # uses the Image module (PIL)
      4 
      5 import matplotlib.pyplot as plt


AttributeError: module 'scipy.misc' has no attribute 'imsave'
1
2
3
4
5
6
7
8
9
10
from scipy import misc
from scipy import ndimage

face = misc.face()
misc.imsave('face.png', face) # First we need to create the PNG file

face = misc.imread('face.png')
print(type(face) )   

print(face.shape, face.dtype)
1
2
3
4
5
face.tofile('face1.raw') # Create raw file
face_from_raw = np.fromfile('face1.raw', dtype=np.uint8)
face_from_raw.shape

face_from_raw.shape = (768, 1024, 3)
1
face_memmap = np.memmap('face.raw', dtype=np.uint8, shape=(768, 1024, 3))
1
2
3
4
5
6
for i in range(10):
    im = np.random.randint(0, 256, 10000).reshape((100, 100))
    misc.imsave('random_%02d.png' % i, im)
from glob import glob
filelist = glob('random*.png')
filelist.sort()
1
2
3
f = misc.face(gray=True)  # retrieve a grayscale image
import matplotlib.pyplot as plt
plt.imshow(f, cmap=plt.cm.gray)  
1
2
3
4
plt.imshow(f, cmap=plt.cm.gray, vmin=30, vmax=200)      

# Remove axes and ticks
plt.axis('off')
1
plt.contour(f, [50, 200])
1
face = misc.face(gray=True)
1
face[0, 40]
1
face[10:13, 20:23]
1
face[100:120] = 255
1
lx, ly = face.shape
1
2
X, Y = np.ogrid[0:lx, 0:ly]
mask = (X - lx / 2) ** 2 + (Y - ly / 2) ** 2 > lx * ly / 4
1
face[mask] = 0
1
face[range(400), range(400)] = 255
1
plt.imshow(face, cmap=plt.cm.gray)  
1

Pandas

creation

  1. Series
    • by passing a list of values, letting pandas create a defualt integer index
    • pd.Series(value,index=list)
  2. DataFrame
    • by passing a NumPy array
    • by passing a dict of objects
    • df.dtypes

viewing data

1
2
3
4
5
6
- df.head();df.tail();df.index;df.columns.
- df.to_numpy():numpy array have one dtype for the entire array, while pandas dataframe have one dtype per column.
- df.describe(): show a quick statistic summary of your data
- df.T : transppsing your data
- df.sort_index(): sorting by an axis
- df.sort_values(): sorting by values

selection

1
2
3
4
5
6
7
8
1. selection by label: df.loc[]; df.at[]
2. selection by position: df.iloc[]; df.iat[]
3. Boolean indexing
    - isin()
4. setting values
    - by indexes
    - by label
    - by position

missing data

1
2
3
4
5
import numpy as np
import pandas as pd

s = pd.Series([1,3,5,np.nan,6,8])
s
1
2
3
4
5
6
7
0    1.0
1    3.0
2    5.0
3    NaN
4    6.0
5    8.0
dtype: float64
1
2
dates = pd.date_range('20130101',periods=6)
dates
1
2
3
DatetimeIndex(['2013-01-01', '2013-01-02', '2013-01-03', '2013-01-04',
               '2013-01-05', '2013-01-06'],
              dtype='datetime64[ns]', freq='D')
1
2
df = pd.DataFrame(np.random.randn(6,4),index=dates,columns=list('ABCD'))
df
A B C D
2013-01-01 0.024195 0.192700 -1.662034 -1.575218
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790
2013-01-04 -1.088395 -0.959165 0.931658 1.973182
2013-01-05 -0.733413 1.216093 0.712507 0.151266
2013-01-06 -0.583591 -0.677252 -0.045332 1.818146 
1
2
3
4
5
6
7
df2 = pd.DataFrame({'A': 1.,
'B': pd.Timestamp('20130102'),
'C': pd.Series(1, index=list(range(4)), dtype='float32'),
 'D': np.array([3] * 4, dtype='int32'),
'E': pd.Categorical(["test", "train", "test", "train"]),
'F': 'foo'})
df2
A B C D E F
0 1.0 2013-01-02 1.0 3 test foo
1 1.0 2013-01-02 1.0 3 train foo
2 1.0 2013-01-02 1.0 3 test foo
3 1.0 2013-01-02 1.0 3 train foo 
1
df2.dtypes
1
2
3
4
5
6
7
A           float64
B    datetime64[ns]
C           float32
D             int32
E          category
F            object
dtype: object
1
df
A B C D
2013-01-01 0.024195 0.192700 -1.662034 -1.575218
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790
2013-01-04 -1.088395 -0.959165 0.931658 1.973182
2013-01-05 -0.733413 1.216093 0.712507 0.151266
2013-01-06 -0.583591 -0.677252 -0.045332 1.818146 
1
df.head()
A B C D
2013-01-01 0.024195 0.192700 -1.662034 -1.575218
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790
2013-01-04 -1.088395 -0.959165 0.931658 1.973182
2013-01-05 -0.733413 1.216093 0.712507 0.151266 
1
df.tail(3)
A B C D
2013-01-04 -1.088395 -0.959165 0.931658 1.973182
2013-01-05 -0.733413 1.216093 0.712507 0.151266
2013-01-06 -0.583591 -0.677252 -0.045332 1.818146 
1
df.index
1
2
3
DatetimeIndex(['2013-01-01', '2013-01-02', '2013-01-03', '2013-01-04',
               '2013-01-05', '2013-01-06'],
              dtype='datetime64[ns]', freq='D')
1
df.columns
1
Index(['A', 'B', 'C', 'D'], dtype='object')
1
df.to_numpy()
1
2
3
4
5
6
array([[ 0.02419536,  0.1927002 , -1.66203397, -1.57521826],
       [ 1.77256918,  2.58352792, -1.61749227,  1.29210627],
       [ 1.40455066,  0.51143439, -0.13094149, -0.18079014],
       [-1.08839452, -0.95916509,  0.93165768,  1.97318197],
       [-0.73341308,  1.21609308,  0.71250658,  0.15126616],
       [-0.58359064, -0.67725157, -0.04533231,  1.81814613]])
1
df2.to_numpy()
1
2
3
4
5
array([[1.0, Timestamp('2013-01-02 00:00:00'), 1.0, 3, 'test', 'foo'],
       [1.0, Timestamp('2013-01-02 00:00:00'), 1.0, 3, 'train', 'foo'],
       [1.0, Timestamp('2013-01-02 00:00:00'), 1.0, 3, 'test', 'foo'],
       [1.0, Timestamp('2013-01-02 00:00:00'), 1.0, 3, 'train', 'foo']],
      dtype=object)
1
df.describe()
A B C D
count 6.000000 6.000000 6.000000 6.000000
mean 0.132653 0.477890 -0.301939 0.579782
std 1.189356 1.300815 1.115929 1.370302
min -1.088395 -0.959165 -1.662034 -1.575218
25% -0.695957 -0.459764 -1.245855 -0.097776
50% -0.279698 0.352067 -0.088137 0.721686
75% 1.059462 1.039928 0.523047 1.686636
max 1.772569 2.583528 0.931658 1.973182 
1
df.T
2013-01-01 2013-01-02 2013-01-03 2013-01-04 2013-01-05 2013-01-06
A 0.024195 1.772569 1.404551 -1.088395 -0.733413 -0.583591
B 0.192700 2.583528 0.511434 -0.959165 1.216093 -0.677252
C -1.662034 -1.617492 -0.130941 0.931658 0.712507 -0.045332
D -1.575218 1.292106 -0.180790 1.973182 0.151266 1.818146 
1
df
A B C D
2013-01-01 0.024195 0.192700 -1.662034 -1.575218
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790
2013-01-04 -1.088395 -0.959165 0.931658 1.973182
2013-01-05 -0.733413 1.216093 0.712507 0.151266
2013-01-06 -0.583591 -0.677252 -0.045332 1.818146 
1
df.sort_index(axis=1,ascending=False)
D C B A
2013-01-01 -1.575218 -1.662034 0.192700 0.024195
2013-01-02 1.292106 -1.617492 2.583528 1.772569
2013-01-03 -0.180790 -0.130941 0.511434 1.404551
2013-01-04 1.973182 0.931658 -0.959165 -1.088395
2013-01-05 0.151266 0.712507 1.216093 -0.733413
2013-01-06 1.818146 -0.045332 -0.677252 -0.583591 
1
df.sort_values(by='B')
A B C D
2013-01-04 -1.088395 -0.959165 0.931658 1.973182
2013-01-06 -0.583591 -0.677252 -0.045332 1.818146
2013-01-01 0.024195 0.192700 -1.662034 -1.575218
2013-01-03 1.404551 0.511434 -0.130941 -0.180790
2013-01-05 -0.733413 1.216093 0.712507 0.151266
2013-01-02 1.772569 2.583528 -1.617492 1.292106 
1
df['A']
1
2
3
4
5
6
7
2013-01-01    0.024195
2013-01-02    1.772569
2013-01-03    1.404551
2013-01-04   -1.088395
2013-01-05   -0.733413
2013-01-06   -0.583591
Freq: D, Name: A, dtype: float64
1
df[0:3]
A B C D
2013-01-01 0.024195 0.192700 -1.662034 -1.575218
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790 
1
df['20130102':'20130104']
A B C D
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790
2013-01-04 -1.088395 -0.959165 0.931658 1.973182 
1
df.loc[dates[0]]
1
2
3
4
5
A    0.024195
B    0.192700
C   -1.662034
D   -1.575218
Name: 2013-01-01 00:00:00, dtype: float64
1
df.loc[:,['A','B']]
A B
2013-01-01 0.024195 0.192700
2013-01-02 1.772569 2.583528
2013-01-03 1.404551 0.511434
2013-01-04 -1.088395 -0.959165
2013-01-05 -0.733413 1.216093
2013-01-06 -0.583591 -0.677252 
1
df.loc['20130102':'20130104',['A','B']]
A B
2013-01-02 1.772569 2.583528
2013-01-03 1.404551 0.511434
2013-01-04 -1.088395 -0.959165 
1
df.loc['20130102', ['A', 'B']]
1
2
3
A    1.772569
B    2.583528
Name: 2013-01-02 00:00:00, dtype: float64
1
df.loc[dates[0], 'A']
1
0.024195360467579575
1
df.at[dates[0], 'A']
1
0.024195360467579575
1
df.iloc[3]
1
2
3
4
5
A   -1.088395
B   -0.959165
C    0.931658
D    1.973182
Name: 2013-01-04 00:00:00, dtype: float64
1
df.iloc[3:5, 0:2]
A B
2013-01-04 -1.088395 -0.959165
2013-01-05 -0.733413 1.216093 
1
df.iloc[[1, 2, 4], [0, 2]]
A C
2013-01-02 1.772569 -1.617492
2013-01-03 1.404551 -0.130941
2013-01-05 -0.733413 0.712507 
1
df.iloc[1:3, :]
A B C D
2013-01-02 1.772569 2.583528 -1.617492 1.292106
2013-01-03 1.404551 0.511434 -0.130941 -0.180790 
1
df.iloc[:, 1:3]
B C
2013-01-01 0.192700 -1.662034
2013-01-02 2.583528 -1.617492
2013-01-03 0.511434 -0.130941
2013-01-04 -0.959165 0.931658
2013-01-05 1.216093 0.712507
2013-01-06 -0.677252 -0.045332

pandas 索引

中括号索引

1
df[0:3]
1
df[0:1]
1
df['A']
1
df.A
1
df['A':'D']
1
df[[0:3],[0:3]]
1

中括号只能索引行与单列

单列的索引，只能是自定义标签，且可以用df.columns_name表示。

行的索引可以用两套索引，但必须带冒号，且顾头不顾尾（按python惯例）。

按标签索引：.loc[[],[]]

同矩阵的索引类似，顾头顾尾，行列的部分的中括号可以酌情省略。

.iat[[],[]],同.iloc[[],[]]，但速度更快。

按位置索引：.iloc[[],[]]

同矩阵的索引类似，顾头不顾尾，行列的部分的中括号可以酌情省略。

.iat[[],[]],同.iloc[[],[]]，但速度更快。

E 布尔索引

1

Scipy

scipy.special

1
2
3
import numpy as np
import matplotlib.pyplot as plt
import scipy as cp
1
from scipy import constants as con
1
con.pi
1
3.141592653589793
1
2
import math
math.pi
1
3.141592653589793
1
con.golden
1
1.618033988749895
1
con.golden_ratio
1
1.618033988749895
1
con.c
1
299792458.0
1
con.speed_of_light
1
299792458.0
1
con.speed_of_sound
1
340.5
1
con.mu_0
1
1.25663706212e-06
1
con.epsilon_0
1
8.8541878128e-12
1
# help(con.physical_constants)
1
con.physical_constants['alpha particle mass']
1
(6.6446573357e-27, 'kg', 2e-36)
1
con.value
1
<function scipy.constants.codata.value(key)>
1
from scipy import special as scs
1
scs.geterr()
1
2
3
4
5
6
7
8
9
{'singular': 'ignore',
 'underflow': 'ignore',
 'overflow': 'ignore',
 'slow': 'ignore',
 'loss': 'ignore',
 'no_result': 'ignore',
 'domain': 'ignore',
 'arg': 'ignore',
 'other': 'ignore'}
1
2
3
4
v = 5
z = np.arange(0,100,0.1)
j = scs.jv(v,z)
plt.plot(z,j)
1
[<matplotlib.lines.Line2D at 0x2b1af1c5280>]

png

1
2
3
4
5
m = 0
v = 10
x = np.linspace(-1,1,100)
a = scs.lpmv(m,v,x)
plt.plot(x,a)
1
[<matplotlib.lines.Line2D at 0x2b1af1d21c0>]

png

1
2


1

scipy.optimize 练习

在数学最优化中，Rosenbrock函数是一个用来测试最优化算法性能的非凸函数，由Howard Harry Rosenbrock在1960年提出。也称为Rosenbrock山谷或Rosenbrock香蕉函数，也简称为香蕉函数。

$f(x,y) = (1-x)^2 +100(y-x^2)^2$

Rosenbrock函数的每个等高线大致呈抛物线形，其全域最小值也位在抛物线形的山谷中（香蕉型山谷）。很容易找到这个山谷，但由于山谷内的值变化不大，要找到全域的最小值相当困难。其全域最小值位于 (x,y)=(1,1)点，数值为f(x,y)=0。有时第二项的系数不同，但不会影响全域最小值的位置。

1
2
def fun_rosenbrock(x):
    return np.array([10 * (x[1] - x[0]**2), (1 - x[0])])
1
2
3
4
5
6
7
8
9
10
import numpy as np
from scipy.optimize import least_squares
x0_rosenbrock = np.array([2, 2])
res_1 = least_squares(fun_rosenbrock, x0_rosenbrock)
print(res_1.x)

print(res_1.cost)

print(res_1.optimality)
print(res_1)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
[1. 1.]
9.866924291084687e-30
8.892886493421953e-14
 active_mask: array([0., 0.])
        cost: 9.866924291084687e-30
         fun: array([4.44089210e-15, 1.11022302e-16])
        grad: array([-8.89288649e-14,  4.44089210e-14])
         jac: array([[-20.00000015,  10.        ],
       [ -1.        ,   0.        ]])
     message: '`gtol` termination condition is satisfied.'
        nfev: 3
        njev: 3
  optimality: 8.892886493421953e-14
      status: 1
     success: True
           x: array([1., 1.])
1
2
3
4
def jac_rosenbrock(x):
    return np.array([
        [-20 * x[0], 10],
        [-1, 0]])
1
2
3
4
5
6
7
8
9
res_2 = least_squares(fun_rosenbrock, x0_rosenbrock, jac_rosenbrock,
                      bounds=([-np.inf, 1.5], np.inf))
print(res_2.x)

print(res_2.cost)

print(res_2.optimality)
print(res_2)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
[1.22437075 1.5       ]
0.025213093946805685
1.5885401433157753e-07
 active_mask: array([ 0, -1])
        cost: 0.025213093946805685
         fun: array([ 0.00916269, -0.22437075])
        grad: array([1.58854014e-07, 9.16268991e-02])
         jac: array([[-24.48741498,  10.        ],
       [ -1.        ,   0.        ]])
     message: '`ftol` termination condition is satisfied.'
        nfev: 9
        njev: 9
  optimality: 1.5885401433157753e-07
      status: 2
     success: True
           x: array([1.22437075, 1.5       ])
1
2
3
4
from scipy.optimize import leastsq
def func(x):
    return 2*(x-3)**2+1
leastsq(func, 0)
1
(array([2.99999999]), 1)
1

python 插值：scipy.interpolate

1
2
from scipy import interpolate as si
import matplotlib.pyplot as plt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 1维插值 interp1d
# 样条插值阶数可选

x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f = si.interp1d(x, y)
f2 = si.interp1d(x, y, kind='cubic')
f3 = si.interp1d(x,y,kind=5)
f4 = si.interp1d(x,y,kind='nearest')
f5 = si.interp1d(x,y,kind='previous')
f6 = si.interp1d(x,y,kind='next')

xnew = np.linspace(0, 10, num=41, endpoint=True)

plt.figure(figsize=(12,9))
plt.subplot(221)
plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')

plt.subplot(222)
plt.plot(x, y, 'o', xnew, f3(xnew), '-', xnew, f4(xnew), '--')
plt.legend(['data', 'Spline5', 'nearest'], loc='best')

plt.subplot(223)
plt.plot(x, y, 'o', xnew, f5(xnew), '-', xnew, f6(xnew), '--')
plt.legend(['data', 'previous', 'next'], loc='best')

plt.tight_layout()
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# 多变量插值：griddata

def func(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2

grid_x, grid_y = np.mgrid[0:1:100j, 0:1:200j]

points = np.random.rand(1000, 2)
values = func(points[:,0], points[:,1])

grid_z0 = si.griddata(points, values, (grid_x, grid_y), method='nearest')
grid_z1 = si.griddata(points, values, (grid_x, grid_y), method='linear')
grid_z2 = si.griddata(points, values, (grid_x, grid_y), method='cubic')

plt.subplot(221)
plt.imshow(func(grid_x, grid_y).T, extent=(0,1,0,1), origin='lower')
plt.plot(points[:,0], points[:,1], 'k.', ms=1)
plt.title('Original')
plt.subplot(222)
plt.imshow(grid_z0.T, extent=(0,1,0,1), origin='lower')
plt.title('Nearest')
plt.subplot(223)
plt.imshow(grid_z1.T, extent=(0,1,0,1), origin='lower')
plt.title('Linear')
plt.subplot(224)
plt.imshow(grid_z2.T, extent=(0,1,0,1), origin='lower')
plt.title('Cubic')
plt.gcf().set_size_inches(6, 6)
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 样条插值
# 1. 计算样条表达式
# 2. 计算插值

x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
y = np.sin(x)
tck = si.splrep(x, y, s=0)
xnew = np.arange(0, 2*np.pi, np.pi/50)
ynew = si.splev(xnew, tck, der=0)

plt.figure()
plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
plt.legend(['Linear', 'Cubic Spline', 'True'])
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('Cubic-spline interpolation')
plt.show()

png

1
2
3
4
5
6
7
yder = si.splev(xnew, tck, der=1)
plt.figure()
plt.plot(xnew, yder, xnew, np.cos(xnew),'--')
plt.legend(['Cubic Spline', 'True'])
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('Derivative estimation from spline')
plt.show()

png

外推法:Scipy.UnivariateSpline

1
2
3
4
5
6
7
8
9
10
11
12
13
from scipy import interpolate
import numpy as np
x1=np.linspace(0,10,20)
y1=np.sin(x1)

sx1=np.linspace(0,12,100)
func1=interpolate.UnivariateSpline(x1,y1,s=0)#强制通过所有点
sy1=func1(sx1)

import matplotlib.pyplot as plt
plt.plot(x1,y1,'o')
plt.plot(sx1,sy1)
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
import numpy as np
from scipy import interpolate
x2=np.linspace(0,20,200)
y2=np.sin(x2)+np.random.normal(loc=0,scale=1,size=len(x2))*0.2

sx2=np.linspace(0,22,2000)
func2=interpolate.UnivariateSpline(x2,y2,s=8)
sy2=func2(sx2)

import matplotlib.pyplot as plt
plt.plot(x2,y2,'.')
plt.plot(sx2,sy2)
plt.show()

png

1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 移动平均算法：一种趋势算法。

def waitui(x,w,n,lr):
    '''
    x:data,list_like;
    w:width
    n:number
    lr:left or right
    '''
    for k in range(n):
        if lr == 'l':
            x0 = w*x[int((w-1)/2)-1] - sum(x[:int(w-1)]) 
            x.insert(0,x0)
        if lr == 'r':
            x0 = w*x[int(0-(w-1)/2)] - sum(x[-int(w-1):])
            x.append(x0)
          
    return x
1
2
3
x2=np.linspace(0,20,20)
y2=2*x2+3 +np.random.normal(loc=0,scale=1,size=len(x2))
plt.plot(x2,y2,'.')
1
[<matplotlib.lines.Line2D at 0x1854963fc40>]

png

1

1
y2
1
2
3
4
array([ 2.13802457,  5.61154911,  6.05997104,  9.16678043,  9.60155538,
       13.1974428 , 15.90229562, 16.02519421, 20.04754014, 22.14220077,
       26.35090094, 26.50786393, 28.57640137, 30.3048343 , 31.92003191,
       33.60209416, 37.30656326, 39.41982125, 41.34787236, 43.59627556])
1
2
3
4
ls = [1,2,3,4,5]
ls = list(y2)
temp = waitui(ls,5,9,'r')
plt.plot(temp)
1
[<matplotlib.lines.Line2D at 0x18546c3ffa0>]

png

1
temp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
[2.1380245690246995,
 5.611549114103028,
 6.059971043082472,
 9.166780428711341,
 9.601555382236675,
 13.197442799881992,
 15.90229562319444,
 16.025194212906964,
 20.04754014250857,
 22.142200770196332,
 26.35090093724794,
 26.50786392690083,
 28.576401370363442,
 30.304834299193807,
 31.920031914422612,
 33.60209415559048,
 37.306563260759276,
 39.41982124700761,
 41.347872357070834,
 43.596275556223034,
 45.84467875537524,
 48.09308195452745,
 50.34148515367964,
 59.33509795028846,
 54.83829155198407,
 97.55795233587588,
 59.33509795028846,
 14.367033967244538,
 347.1307074417702,
 -446.55562185895644,
 1761.3763197085043,
 -3909.0965485533447]
1
plt.plot(temp)
1
[<matplotlib.lines.Line2D at 0x18549450640>]

png

1
3*ls[0] - sum(ls[:2])
1
0
1
ls[:2]
1
[-13, -3]
1
import scipy
1
scipy.interpolate.Akima1DInterpolator.extrapolate??
1
2
3
Type:        member_descriptor
String form: <member 'extrapolate' of '_PPolyBase' objects>
Docstring:   <no docstring>
1
scipy.__version__
1
'1.5.0'
1
2
3
4
5
6
7
8
9
10
import matplotlib.pyplot as plt
from scipy import interpolate
x = np.arange(0, 10)
y = np.exp(-x/3.0)
f = interpolate.interp1d(x, y,kind=5,fill_value='extrapolate')

xnew = np.arange(0, 13, 0.1)
ynew = f(xnew)   # use interpolation function returned by `interp1d`
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.show()

png

光谱平滑 :Scipy.signal.savgol_fliter

  1. 移动平滑
  2. Savitzky-Golay卷积平滑

python 最小二乘相关函数

非线性最小二乘

  • least_squares

变量要有边界

  • 目标函数，非线性最小二乘里面称为残差(测量值与预测值之差)

当我们用一个模型来描述现实中的一系列数据时，模型的预测结果与实际的测量结果总会存在一定偏差，这一偏差就称为残差。非线性最小二乘的目的就是，调整模型的参数，使得总的残差最小。
残差的平方和是最佳的目标函数

  • 损失函数，减少异常值的影响，多种方式可选
  • 代价函数，$0.5\sum \rho(f(x)^2)$
  • 高斯-牛顿法（Gauss-Newton method）
  • 列文伯格-马夸尔特法（Levenberg-Marquardt method）

在给定变量范围内，找到代价函数的局部最小值

参数

  1. 必须给定参数
  • 函数
  • 初值
  1. 其它参数可选
  • 对于非线性模型的最小二乘问题是没有线性解析解的
  • 对没有解析解的问题，通常都是用 逐渐逼近 的思想去解决。而迭代法就正是这种思想的一种重要实现手段。
  • 最小二乘，广义上来说其实是机器学习中的平方损失函数
  • 线性最小二乘有闭式解，可用最小二乘法求解，也可采用迭代法（如梯度下降）求解；非线性最小二乘没有闭式解，只能采用迭代法求解。

线性最小二乘

  • nnls;This is a wrapper for a FORTRAN non-negative least squares solver.
  • lsq_linear;Solve a linear least-squares problem with bounds on the variables.
1
2
3
4
5
from scipy.optimize import nnls

A = np.array([[1, 0], [1, 0], [0, 1]])
b = np.array([2, 1, 1])
nnls(A, b)
1
(array([1.5, 1. ]), 0.7071067811865475)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from scipy.sparse import rand
from scipy.optimize import lsq_linear

np.random.seed(0)

m = 20000
n = 10000

A = rand(m, n, density=1e-4)
b = np.random.randn(m)

lb = np.random.randn(n)
ub = lb + 1

res = lsq_linear(A, b, bounds=(lb, ub), lsmr_tol='auto', verbose=1)
res
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
The relative change of the cost function is less than `tol`.
Number of iterations 15, initial cost 1.5178e+04, final cost 1.1150e+04, first-order optimality 7.15e-08.





 active_mask: array([ 1, -1, -1, ..., -1,  0,  1])
        cost: 11150.33129344828
         fun: array([ 0.00219242, -0.10699277, -0.70574952, ...,  0.76716404,
       -0.40184636,  0.67399679])
     message: 'The relative change of the cost function is less than `tol`.'
         nit: 15
  optimality: 7.146292795202272e-08
      status: 2
     success: True
           x: array([ 0.54987404, -0.05812073,  1.34915899, ...,  1.36754417,
        0.        , -0.96004356])

曲线拟合

  • curve_fit;Use non-linear least squares to fit a function, f, to data.

输入值：

  • 函数：
  • xdata
  • ydata
  • 其他参数可选

返回值：

  • 最优参数；此时，残差平方和最小
  • 参数的估计方差；

To compute one standard deviation errors on the parameters use perr = np.sqrt(np.diag(pcov)).

1
2
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
1
2
def func(x, a, b, c):
    return a * np.exp(-b * x) + c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
np.random.seed(1729)
y_noise = 0.2 * np.random.normal(size=xdata.size)
ydata = y + y_noise
plt.plot(xdata, ydata, 'b-', label='data')

popt, pcov = curve_fit(func, xdata, ydata)
popt

plt.plot(xdata, func(xdata, *popt), 'r-',
         label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))

popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 1., 0.5]))
popt

plt.plot(xdata, func(xdata, *popt), 'g--',
         label='fit: a=%5.3f, b=%5.3f, c=%5.3f' % tuple(popt))

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

png

1
2
xdata = np.linspace(-10, 10, num=20)
ydata = f(xdata) + np.random.randn(xdata.size)
1
plt.plot(xdata,ydata)
1
[<matplotlib.lines.Line2D at 0x2b1b02ed490>]

png

1
2
def f2(x, a, b):
    return a*x**2 + b*np.sin(x)
1
2
3
guess = [2, 2]
params, params_covariance = optimize.curve_fit(f2, xdata, ydata, guess)
params
1
array([1.00278316, 9.10173369])
1
2
a = params[0]
b = params[1]
1
2
plt.plot(xdata,ydata)
plt.plot(x,f2(x,a,b))
1
[<matplotlib.lines.Line2D at 0x2b1b034b700>]

png

Matplotlib

CSDN_matplotlib详解

1
%matplotlib inline

Colormap reference

中文说明见CSDN

记住一个参数：

plt.get_cmap('str')

Reference for colormaps included with Matplotlib.

A reversed version of each of these colormaps is available by appending
_r to the name, e.g., viridis_r.

See :doc:/tutorials/colors/colormaps for an in-depth discussion about
colormaps, including colorblind-friendliness.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import numpy as np
import matplotlib.pyplot as plt


cmaps = [('Perceptually Uniform Sequential', [
            'viridis', 'plasma', 'inferno', 'magma', 'cividis']),
         ('Sequential', [
            'Greys', 'Purples', 'Blues', 'Greens', 'Oranges', 'Reds',
            'YlOrBr', 'YlOrRd', 'OrRd', 'PuRd', 'RdPu', 'BuPu',
            'GnBu', 'PuBu', 'YlGnBu', 'PuBuGn', 'BuGn', 'YlGn']),
         ('Sequential (2)', [
            'binary', 'gist_yarg', 'gist_gray', 'gray', 'bone', 'pink',
            'spring', 'summer', 'autumn', 'winter', 'cool', 'Wistia',
            'hot', 'afmhot', 'gist_heat', 'copper']),
         ('Diverging', [
            'PiYG', 'PRGn', 'BrBG', 'PuOr', 'RdGy', 'RdBu',
            'RdYlBu', 'RdYlGn', 'Spectral', 'coolwarm', 'bwr', 'seismic']),
         ('Cyclic', ['twilight', 'twilight_shifted', 'hsv']),
         ('Qualitative', [
            'Pastel1', 'Pastel2', 'Paired', 'Accent',
            'Dark2', 'Set1', 'Set2', 'Set3',
            'tab10', 'tab20', 'tab20b', 'tab20c']),
         ('Miscellaneous', [
            'flag', 'prism', 'ocean', 'gist_earth', 'terrain', 'gist_stern',
            'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg',
            'gist_rainbow', 'rainbow', 'jet', 'nipy_spectral', 'gist_ncar'])]


gradient = np.linspace(0, 1, 256)
gradient = np.vstack((gradient, gradient))


def plot_color_gradients(cmap_category, cmap_list):
    # Create figure and adjust figure height to number of colormaps
    nrows = len(cmap_list)
    figh = 0.35 + 0.15 + (nrows + (nrows-1)*0.1)*0.22
    fig, axes = plt.subplots(nrows=nrows, figsize=(6.4, figh))
    fig.subplots_adjust(top=1-.35/figh, bottom=.15/figh, left=0.2, right=0.99)

    axes[0].set_title(cmap_category + ' colormaps', fontsize=14)

    for ax, name in zip(axes, cmap_list):
        ax.imshow(gradient, aspect='auto', cmap=plt.get_cmap(name))
        ax.text(-.01, .5, name, va='center', ha='right', fontsize=10,
                transform=ax.transAxes)

    # Turn off *all* ticks & spines, not just the ones with colormaps.
    for ax in axes:
        ax.set_axis_off()


for cmap_category, cmap_list in cmaps:
    plot_color_gradients(cmap_category, cmap_list)

plt.show()

png

png

png

png

png

png

png

References
“”””””””””

The use of the following functions, methods, classes and modules is shown
in this example:

1
2
3
4
5
import matplotlib
matplotlib.colors
matplotlib.axes.Axes.imshow
matplotlib.figure.Figure.text
matplotlib.axes.Axes.set_axis_off
1
<function matplotlib.axes._base._AxesBase.set_axis_off(self)>
1

1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import matplotlib.pyplot as plt
import numpy as np
 
x = np.linspace(-3, 3, 50)
y = 2 * x + 1
plt.figure(num=1, figsize=(8, 5))
plt.plot(x, y)
 
ax = plt.gca()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.spines['bottom'].set_position(('data', 0))
ax.yaxis.set_ticks_position('left')
ax.spines['left'].set_position(('data', 0))
# 添加点
x0 = 1
y0 = 2 * x0 + 1
plt.scatter(x0, y0, s=50, color='b')
plt.plot([x0, x0], [y0, 0], 'k--', lw=2.5)
# 方法一
plt.annotate(r'$2x+1=%s$' % y0, xy=(x0, y0),
             xycoords='data', xytext=(+30, -30),
             textcoords='offset points', fontsize=16,
             arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=.2'))
# 方法二
plt.text(-3.7, 3, r'$This\ is\ the\ some\ text. \mu\ \sigma_i\ \alpha_t$',
         fontdict={'size': 16, 'color': 'r'})
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import matplotlib.pyplot as plt
 
fig = plt.figure()
x = [1, 2, 3, 4, 5, 6, 7]
y = [1, 3, 4, 2, 5, 8, 6]
 
left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
ax1 = fig.add_axes([left, bottom, width, height])
ax1.plot(x, y, 'r')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_title('title')
 
left, bottom, width, height = 0.2, 0.6, 0.25, 0.25
ax2 = fig.add_axes([left, bottom, width, height])
ax2.plot(y, x, 'r')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('title inside 1')
 
plt.axes([.6, .2, .25, .25])
plt.plot(y[::-1], x, 'g')
plt.xlabel('x')
plt.ylabel('y')
plt.title('inside 2')
 
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
import matplotlib.pyplot as plt
import numpy as np
 
x = np.arange(0, 10, 0.1)
y1 = 0.05 * x ** 2
y2 = -1 * y1
 
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(x, y1, 'g-')
ax2.plot(x, y2, 'b-')
 
ax1.set_xlabel('X data')
ax1.set_ylabel('Y1,color=g')
ax2.set_ylabel('Y2,color=b')
 
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import animation
 
fig, ax = plt.subplots()
 
x = np.arange(0, 2 * np.pi, 0.01)
line, = ax.plot(x, np.sin(x))
 
 
def animate(i):
    line.set_ydata(np.sin(x + i / 100))
    return line,
 
 
def init():
    line.set_ydata(np.sin(x))
    return line,
 
 
ani = animation.FuncAnimation(fig=fig, func=animate, frames=100,
                              init_func=init, interval=20, blit=False)
 
plt.show()

png

1

Matplotlib 三维图

1
2
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
ax.plot_surface(X, Y, Z, rstride=8, cstride=8, alpha=0.9, cmap=cm.winter_r)
# 绘制等高线
cset = ax.contourf(X, Y, Z, zdir='z', offset=-100, cmap=cm.coolwarm)
cset = ax.contourf(X, Y, Z, zdir='x', offset=-40, cmap=cm.coolwarm)
cset = ax.contourf(X, Y, Z, zdir='y', offset=40, cmap=cm.coolwarm)

ax.set_xlabel('X')
ax.set_xlim(-40, 40)
ax.set_ylabel('Y')
ax.set_ylim(-40, 40)
ax.set_zlabel('Z')
ax.set_zlim(-100, 100)

plt.show()

png

1

1

1
from scipy import optimize # 优化是寻找最小化或等式的数值解的问题
1
2
3
# 1. 寻找标量函数的最小值
def f(x):
    return x**2 + 10*np.sin(x)
1
2
x = np.arange(-10,10,0.1)
plt.plot(x,f(x)) # 当有%matplotlib inline时，可以省略plt.show()
1
[<matplotlib.lines.Line2D at 0x2b1b28790d0>]

png

1
2
#找到这个函数的最小值的常用有效方式是从给定的初始点开始进行一个梯度下降
optimize.fmin_bfgs(f, 0)
1
2
3
4
5
6
7
8
9
10
11
Optimization terminated successfully.
         Current function value: -7.945823
         Iterations: 5
         Function evaluations: 12
         Gradient evaluations: 6





array([-1.3064401])
1
2
#这个方法的一个可能问题是，如果这个函数有一些局部最低点，算法可能找到这些局部最低点而不是全局最低点，这取决于初始点
optimize.fmin_bfgs(f, 3, disp=0)
1
array([3.83746709])
1
2
3
4
5
6
# 暴力算法找到全局最优点
# 缺点是，当网格变大时，程序会变慢，scipy.optimize.anneal()
# OpenOpt、IPOPT、PyGMO和PyEvolve是关于全局优化的一些有用的包。
grid = (-10, 10, 0.1)
xmin_global = optimize.brute(f, (grid,))
xmin_global
1
array([-1.30641113])
1
2
3
#要找出局部最低点，让我们用scipy.optimize.fminbound将变量限制在(0,10)区间
xmin_local = optimize.fminbound(f, 0, 10)  
xmin_local
1
3.8374671194983834
1

1
# 2. 寻找标量函数的根
1
2
root = optimize.fsolve(f, 1)  # 我们的最初猜想是1
root
1
array([0.])
1
2
root2 = optimize.fsolve(f, -2.5)
root2
1
array([-2.47948183])
1

1
[<matplotlib.lines.Line2D at 0x2b1b034b700>]
1

1

1

直方图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
import numpy as np
import matplotlib.pyplot as plt

# Fixing random state for reproducibility
np.random.seed(19680801)

mu, sigma = 100, 15
x = mu + sigma * np.random.randn(10000)

# the histogram of the data
n, bins, patches = plt.hist(x, np.arange(40,170,0.5), density=True, facecolor='g', alpha=0.75)


plt.xlabel('Smarts')
plt.ylabel('Probability')
# plt.title('Histogram of IQ')
# plt.text(60, .025, r'$\mu=100,\ \sigma=15$')
# plt.xlim(40, 160)
# plt.ylim(0, 0.03)
# plt.grid(True)
plt.show()

png

1
n.shape
1
(259,)
1
len(np.arange(40,170,0.5))
1
260
1
plt.bar(np.arange(40,170,0.5)[:-1],n*0.5)
1
<BarContainer object of 259 artists>

png

1
import pandas as pd
1
2
temp = pd.read_csv(r'C:\Users\hp\Desktop\P59.csv')
temp
OBJECTID pointid grid_code RASTERVALU
0 751834 751834 0.316 0.316
1 751833 751833 0.301 0.301
2 751832 751832 0.277 0.277
3 639012 639012 0.276 0.276
4 639013 639013 0.272 0.272
... ... ... ... ...
1043770 880209 880209 0.001 0.001
1043771 880212 880212 0.001 0.001
1043772 904890 904890 0.001 0.001
1043773 905145 905145 0.001 0.001
1043774 1005453 1005453 0.001 0.001

1043775 rows × 4 columns

1
2
temp = temp[temp.RASTERVALU<0.2]
temp
OBJECTID pointid grid_code RASTERVALU
46 638757 638757 0.199 0.199
47 775890 775890 0.198 0.198
48 170770 170770 0.195 0.195
49 638886 638886 0.194 0.194
50 638756 638756 0.193 0.193
... ... ... ... ...
1043770 880209 880209 0.001 0.001
1043771 880212 880212 0.001 0.001
1043772 904890 904890 0.001 0.001
1043773 905145 905145 0.001 0.001
1043774 1005453 1005453 0.001 0.001

1043729 rows × 4 columns

1
x = temp.RASTERVALU
1
n,bins,patches = plt.hist(x,bins=np.arange(0,0.2,0.001),density=True)

png

1
plt.bar(bins[:-1],n*0.001,width=0.001)
1
<BarContainer object of 199 artists>

png

plt.style

1
plt.style.use('seaborn')
1
print(plt.style.available)
1
['Solarize_Light2', '_classic_test_patch', 'bmh', 'classic', 'dark_background', 'fast', 'fivethirtyeight', 'ggplot', 'grayscale', 'seaborn', 'seaborn-bright', 'seaborn-colorblind', 'seaborn-dark', 'seaborn-dark-palette', 'seaborn-darkgrid', 'seaborn-deep', 'seaborn-muted', 'seaborn-notebook', 'seaborn-paper', 'seaborn-pastel', 'seaborn-poster', 'seaborn-talk', 'seaborn-ticks', 'seaborn-white', 'seaborn-whitegrid', 'tableau-colorblind10']
1

1

1

  • enumerate
  • 生成式
  • 生成器

容器数据类型

  • 字符串
  • 列表

注意多重列表
注意：乘法可重复列表元素

  • 元组：不可修改元素，便于维护，进程安全，创建更快。
  • 集合：和数学上一致，大括号
  • 字典：键值对
  • 构造器语法？
  • 推导式语法？

类和对象
类是对象的蓝图和模板，而对象是类的实例

  • 私有属性和方法：两个下划线开头。
  • 单下划线开头表示属性是受保护的。
1

正则表达式

  1. 速查表：https://www.debuggex.com/cheatsheet/regex/python
  2. 30分钟入门教程：https://deerchao.cn/tutorials/regex/regex.htm
  3. 在线验证网站：https://tool.oschina.net/regex/

正则表示即匹配规则

  • 处理字符串
  • 爬取网页
  • 元字符
  • 限定符
  • 精准匹配
1
2
import re
# - re.findall()
1
2
3
4
5
6
7
# 匹配电话号
str1 = '我的电话号是：13888888888；他的电话号是：13666666666'
result = re.findall('\d\d\d\d\d\d\d\d\d\d\d',str1)
print(result)

result = re.findall('\d{11}',str1)
print(result)
1
2
['13888888888', '13666666666']
['13888888888', '13666666666']
1

1

1

1

1

1

1
2
3
4
5
from scipy.signal import savgol_filter

# savgol_filter(data,n,k)
# 窗口宽度n:为奇数，越大越平滑
# 多项式阶数k:小于n,越小越平滑

用sympy求解方程

示例

1
1-4/(1.5*2.5**2)
1
0.5733333333333333
1
(0.5/2.5)**2+0.05
1
0.09000000000000001
1
0.91*0.43
1
0.3913
1
2
from sympy import *
%matplotlib inline
1
2
3
x, y, z, t = symbols('x y z t')
k, m, n = symbols('k m n', integer=True)
f, g, h = symbols('f g h', cls=Function)
1
2
f = Rational(3,2)*pi + exp(I*x) / (x**2 + y)
f

$\displaystyle \frac{3 \pi}{2} + \frac{e^{i x}}{x^{2} + y}$

1

1
2
x = Symbol('x')
exp(I*x).subs(x,pi).evalf()

$\displaystyle -1.0$

1
2
g = exp(I*x)
g

$\displaystyle e^{i x}$

1
2
h = f + g
h

$\displaystyle e^{i x} + \frac{3 \pi}{2} + \frac{e^{i x}}{x^{2} + y}$

1
exp(I*x).subs(x,pi)

$\displaystyle -1$

1

1

1
expr.args
1
(x, 2*y)
1
exp(pi * sqrt(163)).evalf() # 计算数值解，括号内参数为有效数字位，默认15？

$\displaystyle 2.62537412640769 \cdot 10^{17}$

1
exp(pi * sqrt(163))

$\displaystyle e^{\sqrt{163} \pi}$

1
latex(S('2*4+10',evaluate=False)) # Latex格式
1
'2 \\cdot 4 + 10'
1
latex(exp(x*2)/2)
1
'\\frac{e^{2 x}}{2}'
1
((x+y)**2 * (x+1)).expand() # 展开

$\displaystyle x^{3} + 2 x^{2} y + x^{2} + x y^{2} + 2 x y + y^{2}$

1
((x+y)**2 * (x+1))

$\displaystyle \left(x + 1\right) \left(x + y\right)^{2}$

1
2
a = 1/x + (x*sin(x) - 1)/x
a

$\displaystyle \frac{x \sin{\left(x \right)} - 1}{x} + \frac{1}{x}$

1
simplify(a)

$\displaystyle \sin{\left(x \right)}$

1
solve(Eq(x**3 + 2*x**2 + 4*x + 8, 0), x)
1
[-2, -2*I, 2*I]
1
solve(x**3 + 2*x**2 + 4*x + 8, x)
1
[-2, -2*I, 2*I]
1
solve([Eq(x + 5*y, 2), Eq(-3*x + 6*y, 15)], [x, y])
1
{x: -3, y: 1}
1
solve([x + 5*y - 2, -3*x + 6*y - 15], [x, y])
1
{x: -3, y: 1}
1
2
3
4
y=Function('y')
n=Symbol('n', integer=True)
f=y(n)-2*y(n-1)-5*y(n-2)
f

$\displaystyle y{\left(n \right)} - 5 y{\left(n - 2 \right)} - 2 y{\left(n - 1 \right)}$

1
rsolve(f,y(n),[1,4])

$\displaystyle \left(\frac{1}{2} - \frac{\sqrt{6}}{4}\right) \left(1 - \sqrt{6}\right)^{n} + \left(\frac{1}{2} + \frac{\sqrt{6}}{4}\right) \left(1 + \sqrt{6}\right)^{n}$

1
2
3
a, b = symbols('a b')
s = Sum(6*n**2 + 2**n, (n, a, b))
s

$\displaystyle \sum_{n=a}^{b} \left(2^{n} + 6 n^{2}\right)$

1
s.doit()

$\displaystyle - 2^{a} + 2^{b + 1} - 2 a^{3} + 3 a^{2} - a + 2 b^{3} + 3 b^{2} + b$

1
product(n*(n+1), (n, 1, b))

$\displaystyle {2}^{\left(b\right)} b!$

1
2
3
f=Function('f')
ex=Eq(f(1/x)-3*f(x),x)
ex.subs(x,2)

$\displaystyle f{\left(\frac{1}{2} \right)} - 3 f{\left(2 \right)} = 2$

1
ex.subs(x,Rational(1,2))

$\displaystyle - 3 f{\left(\frac{1}{2} \right)} + f{\left(2 \right)} = \frac{1}{2}$

1
solve([f(Rational(1,2))-3*f(2)-2,-3*f(Rational(1,2))+f(2)-Rational(1,2)])
1
[{f(1/2): -7/16, f(2): -13/16}]
1
limit((sin(x)-x)/x**3, x, 0)

$\displaystyle - \frac{1}{6}$

1
(1/cos(x)).series(x, 0, 6)

$\displaystyle 1 + \frac{x^{2}}{2} + \frac{5 x^{4}}{24} + O\left(x^{6}\right)$

1
diff(cos(x**2)**2 / (1+x), x)

$\displaystyle - \frac{4 x \sin{\left(x^{2} \right)} \cos{\left(x^{2} \right)}}{x + 1} - \frac{\cos^{2}{\left(x^{2} \right)}}{\left(x + 1\right)^{2}}$

1
integrate(x**2 * cos(x), x)

$\displaystyle x^{2} \sin{\left(x \right)} + 2 x \cos{\left(x \right)} - 2 \sin{\left(x \right)}$

1
integrate(x**2 * cos(x), (x, 0, pi/2))

$\displaystyle -2 + \frac{\pi^{2}}{4}$

1
2
f = Function('f')
dsolve(Eq(Derivative(f(x),x,x) + 9*f(x), 1), f(x))

$\displaystyle f{\left(x \right)} = C_{1} \sin{\left(3 x \right)} + C_{2} \cos{\left(3 x \right)} + \frac{1}{9}$

1
2
f = Function("f")
Eq(f(x).diff(x, x) + 9*f(x), 1)

$\displaystyle 9 f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = 1$

1
dsolve(_, f(x))

$\displaystyle f{\left(x \right)} = C_{1} \sin{\left(3 x \right)} + C_{2} \cos{\left(3 x \right)} + \frac{1}{9}$

1
2
3
4
5
6
7
8
9
def planck(w, T):
    h = 6.626e-34
    c = 3.0e+8
    k = 1.38e-23

    a = 2.0*h*c**2
    b = h*c/(w*k*T)
    intensity = a/ ( (w**5) * (np.exp(b) - 1.0) )
    return intensity
1
2
3
4
5
6
7
8
9
h = 6.626e-34
c = 3.0e+8
k = 1.38e-23
T = symbols('T')
lamda = symbols('lamda')
k, h, c = symbols('k h c', integer=True)
P = symbols('P', cls=Function)
P = 2*h*c**2/((lamda**5) * (exp(h*c/(lamda*k*T)) - 1))
P

$\displaystyle \frac{2 c^{2} h}{\lambda^{5} \left(e^{\frac{c h}{T k \lambda}} - 1\right)}$

1
2
result = P.subs(c,3.0e+8).subs(h,6.626e-34).subs(k,1.38e-23).subs(lamda,500e-9).subs(T,5777)
result.evalf(5)

$\displaystyle 2.6237 \cdot 10^{13}$

1
# solve(Eq(P.subs(c,3.0e+8).subs(h,6.626e-34).subs(k,1.38e-23).subs(lamda,500e-9), 2.6237e13), T)
1
2
3
4
t = symbols('t',cls=Function)
Ra = symbols('I')
t = h*c/(lamda*k*ln((2*h*c**2/(Ra*lamda**5))+1))
t

$\displaystyle \frac{c h}{k \lambda \log{\left(1 + \frac{2 c^{2} h}{I \lambda^{5}} \right)}}$

1
t.subs(c,3.0e+8).subs(h,6.626e-34).subs(k,1.38e-23).subs(lamda,500e-9).subs(Ra,2.6237e13)

$\displaystyle 5777.00277916166$

1
2
3
4
r0 = symbols('r0',cls=Function)
w = symbols('w')
r0 = (1-sqrt(1-w))/(1+sqrt(1-w))
r0

$\displaystyle \frac{1 - \sqrt{1 - w}}{\sqrt{1 - w} + 1}$

1
solve(Eq(r0,0),w)
1
[0]
1
2
3
4
P = symbols('P',cls=Function)
b,c,g = symbols('b c g')
P = 1+b*cos(pi*g/180)+c*(1.5*cos(pi*g/180)**2-0.5)
P

$\displaystyle b \cos{\left(\frac{\pi g}{180} \right)} + c \left(1.5 \cos^{2}{\left(\frac{\pi g}{180} \right)} - 0.5\right) + 1$

1
cos(x).subs(x,pi).evalf()

$\displaystyle -1.0$

1
solve(Eq(P.subs(b,-0.17).subs(c,0.7),1),g)
1
[48.3981512828419, 120.135426949009, 239.864573050991, 311.601848717158]
1
2
def r1(w):
    return (1-np.sqrt(1-w))/(1+np.sqrt(1-w))
1
2
3
f = symbols('f',cls=Function)
f = r0 + P
f
1
2
3
4
5
6
7
8
9
10
11
---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

<ipython-input-76-b0e877824aad> in <module>
      1 f = symbols('f',cls=Function)
----> 2 f = r0 + P
      3 f


TypeError: unsupported operand type(s) for +: 'function' and 'Add'
1
import sympy as sy
1
2
3
4
5
6
7
8
9
10
# 1.表达式与表达式求值

# 多项式求解:
# 定义变量：
x = sy.Symbol('x')
fx = 5*x+4
print(fx)
# 使用evalf函数传值
y1 = fx.evalf(subs={x:6})
print(y1)
1
2
5*x + 4
34.0000000000000
1
2
3
4
5
6
#多元表达式
x=sy.Symbol('x')
y=sy.Symbol('y')
fx=x*x+y*y
result=fx.evalf(subs={x:3,y:4})
print(result)
1
25.0000000000000
1
2
3
4
5
6
7
8
9
10
# 2.函数方程求解

#解方程 有限解
#定义变量
x=sy.Symbol('x')
y=sy.Symbol('y')
fx=x*3+9
#可求解直接给出解向量
print(sy.solve(fx,x))
print(sy.solve(Eq(fx,1),x))
1
2
[-3]
[-8/3]
1
2
3
4
5
6
7
#解方程 无穷多解
#定义变量
x=sy.Symbol('x')
y=sy.Symbol('y')
fx=x*3+y**2
#得到是x与y的关系式，
print(sy.solve(fx,x,y))
1
[(-y**2/3, y)]
1
2
3
4
5
6
7
8
#解方程组
#定义变量
x=sy.Symbol('x')
y=sy.Symbol('y')
f1=x+y-3
f2=x-y+5
print(sy.solve([f1,f2],[x,y]))
# print(sy.solve(Eqs([f1,f2],[0,0]),[x,y]))
1
{x: -1, y: 4}
1
2
3
4
5
6
7
# 求和
#定义变量
n=sy.Symbol('n')
f=2*n
#前面参数放函数，后面放变量的变化范围
s=sy.summation(f,(n,1,100))
print(s)
1
10100
1
2
3
4
5
6
7
8
9
# 解带有求和式的方程

#解释一下，i可以看做是循环变量，就是x自己加五次
#先定义变量，再写出方程
x=sy.Symbol('x')
i=sy.Symbol('i')
f=sy.summation(x,(i,1,5))+10*x-15
result=sy.solve(f,x)
print(result)
1
[1]
1
2
3
4
5
6
7
8
9
10
11
12
# 求极限
#求极限使用limit方法
#定义变量与函数
x=sy.Symbol('x')
f1=sy.sin(x)/x
f2=(1+x)**(1/x)
f3=(1+1/x)**x
#三个参数是 函数，变量，趋向值
lim1=sy.limit(f1,x,0)
lim2=sy.limit(f2,x,0)
lim3=sy.limit(f3,x,sy.oo) # 趋于无穷大
print(lim1,lim2,lim3)
1
1 E E
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
# 求导

#求导使用diff方法
x=sy.Symbol('x')
f1=2*x**4+3*x+6
#参数是函数与变量
f1_=sy.diff(f,x)
print(f1_)
 
f2=sy.sin(x)
f2_=sy.diff(f2,x)
print(f2_)
 
#求偏导
y=sy.Symbol('y')
f3=2*x**2+3*y**4+2*y
#对x，y分别求导，即偏导
f3_x=sy.diff(f3,x)
f3_y=sy.diff(f3,y)
print(f3_x)
print(f3_y)
1
2
3
4
15
cos(x)
4*x
12*y**3 + 2
1
2
3
4
5
6
7
8
# 求定积分

#求定积分用 integrate方法
x=sy.Symbol('x')
f=2*x
#参数传入 函数，积分变量和范围
result=sy.integrate(f,(x,0,1))
print(result)
1
1
1
2
3
4
5
6
7
# 求定积分

from scipy import integrate
def f(x):
    return x + 1
v, err = integrate.quad(f, 1, 2)# err为误差
print(v)
1
2.5
1
2
3
4
5
6
7
#求多重积分，
# 先求里面的积分，再求外面的
x,t=sy.symbols('x t')
f1=2*t
f2=sy.integrate(f1,(t,0,x))
result=sy.integrate(f2,(x,0,3))
print(result)
1
9
1

1

可视化图表交互—Bokeh

Bokeh 是 Python 的交互式可视图库，用于生成在浏览器里显 示的大规模数据集高性能可视图。

github有翻译的教程

1
2
3
4
5
6
7
8
import numpy as np
import pandas as pd
import os
import os.path as op
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import bokeh
  1. 介绍与设置
  • Bokeh是什么？Bokeh是用于现代浏览器端展示的交互式可视化库。
  • 我能用Bokeh 做什么？
  • 进行环境设置
    • 依赖包
    • 示例数据：bokeh sampledata
1
2
3
4
5
from bokeh.io import output_file,show,output_notebook
from bokeh.plotting import figure
from bokeh.models import HoverTool
# from bokeh.charts import Scatter,Bar,BoxPlot,Chord #No module named 'bokeh.charts'
# from bokeh.layouts import raw
1

1
output_notebook()
Loading BokehJS ... 
1
2
3
4
5
6
from numpy import cos, linspace
x = linspace(-6, 6, 100)
y = cos(x)
p = figure(width=500, height=500)
p.circle(x, y) #, size=7, color="firebrick", alpha=0.5
show(p)
1

1
2
3
4
5
6
7
8
9
from bokeh.sampledata.autompg import autompg

grouped = autompg.groupby("yr")
mpg = grouped["mpg"]
avg = mpg.mean()
std = mpg.std()
years = list(grouped.groups.keys())
american = autompg[autompg["origin"]==1]
japanese = autompg[autompg["origin"]==3]
1
autompg
mpg cyl displ hp weight accel yr origin name
0 18.0 8 307.0 130 3504 12.0 70 1 chevrolet chevelle malibu
1 15.0 8 350.0 165 3693 11.5 70 1 buick skylark 320
2 18.0 8 318.0 150 3436 11.0 70 1 plymouth satellite
3 16.0 8 304.0 150 3433 12.0 70 1 amc rebel sst
4 17.0 8 302.0 140 3449 10.5 70 1 ford torino
... ... ... ... ... ... ... ... ... ...
387 27.0 4 140.0 86 2790 15.6 82 1 ford mustang gl
388 44.0 4 97.0 52 2130 24.6 82 2 vw pickup
389 32.0 4 135.0 84 2295 11.6 82 1 dodge rampage
390 28.0 4 120.0 79 2625 18.6 82 1 ford ranger
391 31.0 4 119.0 82 2720 19.4 82 1 chevy s-10

392 rows × 9 columns

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
p = figure(title="MPG by Year (Japan and US)",tools="")

p.vbar(x=years, bottom=avg-std, top=avg+std, width=0.8, 
       fill_alpha=0.2, line_color=None, legend_label="MPG 1 stddev")

p.circle(x=japanese["yr"], y=japanese["mpg"], size=10, alpha=0.5,
         color="red", legend_label="Japanese")

p.triangle(x=american["yr"], y=american["mpg"], size=10, alpha=0.3,
           color="blue", legend_label="American")

p.legend.location = "top_left"

p.add_tools(HoverTool(tooltips=[("MPG", "@mpg_mean"), ("Cyl, Mfr", "@cyl_mfr")]))
show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
from bokeh.models import ColumnDataSource
from bokeh.layouts import gridplot

source = ColumnDataSource(autompg)

options = dict(plot_width=300, plot_height=300,
               tools="pan,wheel_zoom,box_zoom,box_select,lasso_select")

p1 = figure(title="MPG by Year", **options)
p1.circle("yr", "mpg", color="blue", source=source)

p2 = figure(title="HP vs. Displacement", **options)
p2.circle("hp", "displ", color="green", source=source)

p3 = figure(title="MPG vs. Displacement", **options)
p3.circle("mpg", "displ", size="cyl", line_color="red", fill_color=None, source=source)

p = gridplot([[ p1, p2, p3]], toolbar_location="right")

show(p)
1
2
3
4
5
6
# Plot a complex chart with intearctive hover in a few lines of code

from bokeh.models import ColumnDataSource, HoverTool
from bokeh.plotting import figure
from bokeh.sampledata.autompg import autompg_clean as df
from bokeh.transform import factor_cmap
1
2
df.cyl = df.cyl.astype(str)
df.yr = df.yr.astype(str)
1
2
group = df.groupby(['cyl', 'mfr'])
source = ColumnDataSource(group)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
p = figure(plot_width=800, plot_height=300, title="Mean MPG by # Cylinders and Manufacturer",
           x_range=group, toolbar_location=None, tools="")

p.xgrid.grid_line_color = None
p.xaxis.axis_label = "Manufacturer grouped by # Cylinders"
p.xaxis.major_label_orientation = 1.2

index_cmap = factor_cmap('cyl_mfr', palette=['#2b83ba', '#abdda4', '#ffffbf', '#fdae61', '#d7191c'], 
                         factors=sorted(df.cyl.unique()), end=1)

p.vbar(x='cyl_mfr', top='mpg_mean', width=1, source=source,
       line_color="white", fill_color=index_cmap, 
       hover_line_color="black", hover_fill_color=index_cmap)

p.add_tools(HoverTool(tooltips=[("MPG", "@mpg_mean"), ("Cyl, Mfr", "@cyl_mfr")]))

show(p)
  1. 基本绘图

  • 导入与设置：

    • figure
    • output_file/output_notebook/output_server
    • show/save
    1
    2
    from bokeh.io import output_notebook, show
    from bokeh.plotting import figure # mid-level接口：从简单的默认图表开始（使用合理的默认工具、网格和轴），添加标记和其他视觉属性直接与数据绑定的图形。
  • 基本的散点图：
1
p = figure(plot_width=400, plot_height=400)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
- p.circle
- p.square
- asterisk()
- circle()
- circle_cross()
- circle_x()
- cross()
- diamond()
- diamond_cross()
- inverted_triangle()
- square()
- square_cross()
- square_x()
- triangle()
- x()

  • 基本的线形图

    • p.line()

  • 面图

    • p.image_rgba()
1
2
3
4
5
6
7
# create a new plot with default tools, using figure
p = figure(plot_width=400, plot_height=400)

# add a circle renderer with a size, color, and alpha
p.circle([1, 2, 3, 4, 5], [6, 7, 2, 4, 5], size=15, line_color="navy", fill_color="orange", fill_alpha=0.5)

show(p) # show the results
1
2
3
4
5
6
7
# create a new plot using figure
p = figure(plot_width=400, plot_height=400)

# add a square renderer with a size, color, alpha, and sizes
p.square([1, 2, 3, 4, 5], [6, 7, 2, 4, 5], size=[10, 15, 20, 25, 30], color="firebrick", alpha=0.6)

show(p) # show the results
1
2
3
4
5
6
7
# create a new plot (with a title) using figure
p = figure(plot_width=400, plot_height=400, title="My Line Plot")

# add a line renderer
p.line([1, 2, 3, 4, 5], [6, 7, 2, 4, 5], line_width=2)

show(p) # show the results
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from __future__ import division
import numpy as np

# set up some data
N = 20
img = np.empty((N,N), dtype=np.uint32)
view = img.view(dtype=np.uint8).reshape((N, N, 4))
for i in range(N):
    for j in range(N):
        view[i, j, 0] = int(i/N*255) # red
        view[i, j, 1] = 158          # green
        view[i, j, 2] = int(j/N*255) # blue
        view[i, j, 3] = 255          # alpha
      
# create a new plot (with a fixed range) using figure
p = figure(x_range=[0,10], y_range=[0,10])

# add an RGBA image renderer
p.image_rgba(image=[img], x=[0], y=[0], dw=[10], dh=[10])

show(p) # show the results
1
2
3
4
5
6
7
8
9
10
11
12
# set up some data
x = [1, 2, 3, 4, 5]
y = [6, 7, 8, 7, 3]

# create a new plot with figure
p = figure(plot_width=400, plot_height=400)

# add both a line and circles on the same plot
p.line(x, y, line_width=2)
p.circle(x, y, fill_color="white", size=8)

show(p) # show the results
  1. 样式与主题

  • 颜色与属性

    • 颜色：
      • 颜色名称：‘green’等
      • 16进制RGB（A）值
      • 整数（r,g,b）(0-255)
      • (r,g,b,a)(0-255;0-1)
    • 属性：
      • 三个API:models，plotting,charts
      • 三种视觉样式属性：
        • line:line_color,line_alpha,line_width,line_dash
        • fill:fill_color,fill_alpha
        • text:text_font,text_font_size,text_color,text_alpha

  • Plots(框图)

    • outline
      • outline_line_width
      • outline_line_alpha
      • outline_line_color
    • border

  • Glyphs(标记符号)

    • 选择和非选择 视觉效果
  • Axes(轴)
  • Grid(网格)
1
2
3
4
5
6
7
8
9
# create a new plot with a title
p = figure(plot_width=400, plot_height=400)
p.outline_line_width = 7
p.outline_line_alpha = 0.3
p.outline_line_color = "navy"

p.circle([1,2,3,4,5], [2,5,8,2,7], size=10)

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
p = figure(plot_width=400, plot_height=400)

# keep a reference to the returned GlyphRenderer
r = p.circle([1,2,3,4,5], [2,5,8,2,7])

r.glyph.size = 50
r.glyph.fill_alpha = 0.2
r.glyph.line_color = "firebrick"
r.glyph.line_dash = [5, 1]
r.glyph.line_width = 2

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
p = figure(plot_width=400, plot_height=400, tools="tap", title="Select a circle")
renderer = p.circle([1, 2, 3, 4, 5], [2, 5, 8, 2, 7], size=50,

                    # set visual properties for selected glyphs
                    selection_color="firebrick",

                    # set visual properties for non-selected glyphs
                    nonselection_fill_alpha=0.2,
                    nonselection_fill_color="grey",
                    nonselection_line_color="firebrick",
                    nonselection_line_alpha=1.0)

show(p)
  1. 添加注释
  • Span :Spans是“无限”的垂直线或水平线。
  • Box Annotations（Box注释）
  • Label（标签）
  • LabelSet（标签集）
  • Arrows（箭头）
  • Legends（图例）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np
from bokeh.models.annotations import Span

x = np.linspace(0, 20, 200)
y = np.sin(x)

p = figure(y_range=(-2, 2))
p.line(x, y)

upper = Span(location=1, dimension='width', line_color='olive', line_width=4)
p.add_layout(upper)

lower = Span(location=-1, dimension='width', line_color='firebrick', line_width=4)
p.add_layout(lower)

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
import numpy as np
from bokeh.models.annotations import BoxAnnotation

x = np.linspace(0, 20, 200)
y = np.sin(x)

p = figure(y_range=(-2, 2))
p.line(x, y)

# region that always fills the top of the plot
upper = BoxAnnotation(bottom=1, fill_alpha=0.1, fill_color='olive')
p.add_layout(upper)

# region that always fills the bottom of the plot
lower = BoxAnnotation(top=-1, fill_alpha=0.1, fill_color='firebrick')
p.add_layout(lower)

# a finite region
center = BoxAnnotation(top=0.6, bottom=-0.3, left=7, right=12, fill_alpha=0.1, fill_color='navy')
p.add_layout(center)

show(p)
1
2
3
4
5
6
7
8
9
10
from bokeh.models.annotations import Label
from bokeh.plotting import figure

p = figure(x_range=(0,10), y_range=(0,10))
p.circle([2, 5, 8], [4, 7, 6], color="olive", size=10)

label = Label(x=5, y=7, x_offset=12, text="Second Point", text_baseline="middle")
p.add_layout(label)

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from bokeh.plotting import figure
from bokeh.models import ColumnDataSource, LabelSet


source = ColumnDataSource(data=dict(
    temp=[166, 171, 172, 168, 174, 162],
    pressure=[165, 189, 220, 141, 260, 174],
    names=['A', 'B', 'C', 'D', 'E', 'F']))

p = figure(x_range=(160, 175))
p.scatter(x='temp', y='pressure', size=8, source=source)
p.xaxis.axis_label = 'Temperature (C)'
p.yaxis.axis_label = 'Pressure (lbs)'

labels = LabelSet(x='temp', y='pressure', text='names', level='glyph',
                  x_offset=5, y_offset=5, source=source, render_mode='canvas')


p.add_layout(labels)

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from bokeh.models.annotations import Arrow
from bokeh.models.arrow_heads import OpenHead, NormalHead, VeeHead

p = figure(plot_width=600, plot_height=600)

p.circle(x=[0, 1, 0.5], y=[0, 0, 0.7], radius=0.1,
         color=["navy", "yellow", "red"], fill_alpha=0.1)

p.add_layout(Arrow(end=OpenHead(line_color="firebrick", line_width=4),
                   x_start=0, y_start=0, x_end=1, y_end=0))

p.add_layout(Arrow(end=NormalHead(fill_color="orange"),
                   x_start=1, y_start=0, x_end=0.5, y_end=0.7))

p.add_layout(Arrow(end=VeeHead(size=35), line_color="red",
                   x_start=0.5, y_start=0.7, x_end=0, y_end=0))

show(p)
1
2
3
4
5
6
7
8
9
10
11
import numpy as np

x = np.linspace(0, 4*np.pi, 100)
y = np.sin(x)

p = figure(height=400)

p.circle(x, y, legend_label="sin(x)")
p.line(x, 2*y, legend_label="2*sin(x)", line_dash=[4, 4], line_color="orange", line_width=2)

show(p)
1
2
p.circle(x, y, legend_label="sin(x)")
p.line(x, y, legend_label="sin(x)", line_dash=[4, 4], line_color="orange", line_width=2)
GlyphRenderer(
id = '5269', …)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
from bokeh.sampledata.autompg import autompg
from bokeh.models import LinearColorMapper, ColorBar
from bokeh.palettes import Viridis256

source = ColumnDataSource(autompg)
color_mapper = LinearColorMapper(palette=Viridis256, low=autompg.weight.min(), high=autompg.weight.max())

p = figure(x_axis_label='Year', y_axis_label='MPG', tools='', toolbar_location='above')
p.circle(x='yr', y='mpg', color={'field': 'weight', 'transform': color_mapper}, size=20, alpha=0.6, source=source)

color_bar = ColorBar(color_mapper=color_mapper, label_standoff=12, location=(0,0), title='Weight')
p.add_layout(color_bar, 'right')

show(p)
  1. 数据源和数据转换
  • Bokeh的ColumnDataSource类型
    • 通过 Python Dicts 创建
    • 通过 Pandas DataFrames 创建
1
2
3
4
5
6
7
8
from bokeh.models import ColumnDataSource
source = ColumnDataSource(data={
    'x' : [1, 2, 3, 4, 5],
    'y' : [3, 7, 8, 5, 1],
})
p = figure(plot_width=400, plot_height=400)
p.circle('x', 'y', size=20, source=source)
show(p)
1
2
3
4
5
6
from bokeh.sampledata.iris import flowers as df

source = ColumnDataSource(df)
p = figure(plot_width=400, plot_height=400)
p.circle('petal_length', 'petal_width', source=source)
show(p)
  1. 演示布局：把多个图表（plot）放在一个布局（layout）里
  • 行与列：bokeh.layouts 模块提供了 row 和 column 功能以便在垂直和水平布局中排列图表
  • 网格图：Bokeh还在 bokeh.layouts 模块中提供了一个 gridplot 布局，以便把图表排列成一个网格图
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from bokeh.layouts import row

x = list(range(11))
y0, y1, y2 = x, [10-i for i in x], [abs(i-5) for i in x]

# create a new plot
s1 = figure(width=250, plot_height=250)
s1.circle(x, y0, size=10, color="navy", alpha=0.5)

# create another one
s2 = figure(width=250, height=250)
s2.triangle(x, y1, size=10, color="firebrick", alpha=0.5)

# create and another
s3 = figure(width=250, height=250)
s3.square(x, y2, size=10, color="olive", alpha=0.5)

# show the results in a row
show(row(s1, s2, s3))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from bokeh.layouts import gridplot

# create a new plot
s1 = figure(width=250, plot_height=250)
s1.circle(x, y0, size=10, color="navy", alpha=0.5)

# create another one
s2 = figure(width=250, height=250)
s2.triangle(x, y1, size=10, color="firebrick", alpha=0.5)

# create and another
s3 = figure(width=250, height=250)
s3.square(x, y2, size=10, color="olive", alpha=0.5)

# put all the plots in a gridplot
p = gridplot([[s1, s2], [s3, None]], toolbar_location=None)

# show the results
show(p)
  1. 关联和交互：把不同的图表联系起来，把图表和控件（widget）联系起来
  • Linked Interactions：
    • 不同的Bokeh图表之间是可以有联系的，两个（或更多）图的范围可能是由联系的，当一个图被平移（或放大，或其范围变化），其它的图也相应更新保持一致。还可以将两个图的选择（selection）联系起来，当在一个图上选择某项时，另一个图中对应的项也被选中。
  • Linked panning（相连的平移）
  • Linked brushing
    • 相连的选择（selection）通过类似的方式完成，在图块之间共享数据源。注意，通常 bokeh.plotting 和 bokeh.charts 会自动创建一个缺省数据源。然而，为了共享数据源，我们必须手动创建它们并显式地传递。
  • Hover Tools（悬停工具）
    • Bokeh有一个悬停的工具，当用户将鼠标悬停在一个特定的标记符号（glyph）上时，可以在一个弹出框里显示更多的信息。基本的悬停工具配置相当于提供一个 (name, format) 元组列表。
  • Widgets（控件）
    • Bokeh直接集成了一个小的基本控件集。这些控件和Bokeh服务器或 CustomJS 模型协作使用可以加入更多的互动的功能。您可以在用户指南的 Adding Widgets 部分看到完整的widget列表和示例代码。
    • CustomJS Callbacks（CustomJS回调）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from bokeh.layouts import gridplot

x = list(range(11))
y0, y1, y2 = x, [10-i for i in x], [abs(i-5) for i in x]

plot_options = dict(width=250, plot_height=250, tools='pan,wheel_zoom')

# create a new plot
s1 = figure(**plot_options)
s1.circle(x, y0, size=10, color="navy")

# create a new plot and share both ranges
s2 = figure(x_range=s1.x_range, y_range=s1.y_range, **plot_options)
s2.triangle(x, y1, size=10, color="firebrick")

# create a new plot and share only one range
s3 = figure(x_range=s1.x_range, **plot_options)
s3.square(x, y2, size=10, color="olive")

p = gridplot([[s1, s2, s3]])

# show the results
show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from bokeh.models import ColumnDataSource

x = list(range(-20, 21))
y0, y1 = [abs(xx) for xx in x], [xx**2 for xx in x]

# create a column data source for the plots to share
source = ColumnDataSource(data=dict(x=x, y0=y0, y1=y1))

TOOLS = "box_select,lasso_select,help"

# create a new plot and add a renderer
left = figure(tools=TOOLS, width=300, height=300)
left.circle('x', 'y0', source=source)

# create another new plot and add a renderer
right = figure(tools=TOOLS, width=300, height=300)
right.circle('x', 'y1', source=source)

p = gridplot([[left, right]])

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
from bokeh.models import HoverTool

source = ColumnDataSource(
        data=dict(
            x=[1, 2, 3, 4, 5],
            y=[2, 5, 8, 2, 7],
            desc=['A', 'b', 'C', 'd', 'E'],
        )
    )

hover = HoverTool(
        tooltips=[
            ("index", "$index"),
            ("(x,y)", "($x, $y)"),
            ("desc", "@desc"),
        ]
    )

p = figure(plot_width=300, plot_height=300, tools=[hover], title="Mouse over the dots")

p.circle('x', 'y', size=20, source=source)

show(p)
1
2
3
4
5
6
7
from bokeh.layouts import widgetbox
from bokeh.models.widgets import Slider


slider = Slider(start=0, end=10, value=1, step=.1, title="foo")

show(widgetbox(slider))
1
BokehDeprecationWarning: 'WidgetBox' is deprecated and will be removed in Bokeh 3.0, use 'bokeh.models.Column' instead
1
2
3
4
5
6
7
8
9
10
from bokeh.models import TapTool, CustomJS, ColumnDataSource

callback = CustomJS(code="alert('hello world')")
tap = TapTool(callback=callback)

p = figure(plot_width=600, plot_height=300, tools=[tap])

p.circle(x=[1, 2, 3, 4, 5], y=[2, 5, 8, 2, 7], size=20)

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
from bokeh.layouts import column
from bokeh.models import CustomJS, ColumnDataSource, Slider

x = [x*0.005 for x in range(0, 201)]

source = ColumnDataSource(data=dict(x=x, y=x))

plot = figure(plot_width=400, plot_height=400)
plot.line('x', 'y', source=source, line_width=3, line_alpha=0.6)

slider = Slider(start=0.1, end=6, value=1, step=.1, title="power")

update_curve = CustomJS(args=dict(source=source, slider=slider), code="""
    var data = source.get('data');
    var f = slider.value;
    x = data['x']
    y = data['y']
    for (i = 0; i < x.length; i++) {
        y[i] = Math.pow(x[i], f)
    }
    source.change.emit();
""")
slider.js_on_change('value', update_curve)


show(column(slider, plot))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
from random import random

x = [random() for x in range(500)]
y = [random() for y in range(500)]
color = ["navy"] * len(x)

s = ColumnDataSource(data=dict(x=x, y=y, color=color))
p = figure(plot_width=400, plot_height=400, tools="lasso_select", title="Select Here")
p.circle('x', 'y', color='color', size=8, source=s, alpha=0.4)

s2 = ColumnDataSource(data=dict(xm=[0,1],ym=[0.5, 0.5]))
p.line(x='xm', y='ym', color="orange", line_width=5, alpha=0.6, source=s2)

s.callback = CustomJS(args=dict(s2=s2), code="""
    var inds = cb_obj.get('selected')['1d'].indices;
    var d = cb_obj.get('data');
    var ym = 0
  
    if (inds.length == 0) { return; }
  
    for (i = 0; i < d['color'].length; i++) {
        d['color'][i] = "navy"
    }
    for (i = 0; i < inds.length; i++) {
        d['color'][inds[i]] = "firebrick"
        ym += d['y'][inds[i]]
    }
  
    ym /= inds.length
    s2.get('data')['ym'] = [ym, ym]
  
    cb_obj.trigger('change');
    s2.trigger('change');
""")

show(p)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
---------------------------------------------------------------------------

AttributeError                            Traceback (most recent call last)

<ipython-input-47-0b628e0d1f9d> in <module>
     12 p.line(x='xm', y='ym', color="orange", line_width=5, alpha=0.6, source=s2)
     13 
---> 14 s.callback = CustomJS(args=dict(s2=s2), code="""
     15     var inds = cb_obj.get('selected')['1d'].indices;
     16     var d = cb_obj.get('data');


E:\ProgramFiles\Anaconda\lib\site-packages\bokeh\core\has_props.py in __setattr__(self, name, value)
    283                 matches, text = props, "possible"
    284 
--> 285             raise AttributeError("unexpected attribute '%s' to %s, %s attributes are %s" %
    286                 (name, self.__class__.__name__, text, nice_join(matches)))
    287 


AttributeError: unexpected attribute 'callback' to ColumnDataSource, similar attributes are js_event_callbacks
1
$ bokeh sampledata
1
from bokeh.plotting import figure
1
2
3
4
5
6
7
8
import pandas as pd

from bokeh.palettes import Spectral4
from bokeh.plotting import figure, output_file, show
from bokeh.sampledata.stocks import AAPL, GOOG, IBM, MSFT

p = figure(plot_width=800, plot_height=250, x_axis_type="datetime")
p.title.text = 'Click on legend entries to hide the corresponding lines'
1
2
3
4
5
6
7
8
9
10
11
for data, name, color in zip([AAPL, IBM, MSFT, GOOG], ["AAPL", "IBM", "MSFT", "GOOG"], Spectral4):
    df = pd.DataFrame(data) # 将字典转换为DataFrame
    df['date'] = pd.to_datetime(df['date']) 
    p.line(df['date'], df['close'], line_width=2, color=color, alpha=0.8, legend_label=name)

p.legend.location = "top_left"
p.legend.click_policy="hide"

# output_file("interactive_legend.html", title="interactive_legend.py example")

show(p)
1
2
3
4
5
6
7
from bokeh.io import show
from bokeh.models import Button, CustomJS

button = Button(label="Foo", button_type="success")
button.js_on_click(CustomJS(code="console.log('button: click!', this.toString())"))

show(button)
1
2
3
4
5
6
7
8
9
10
11
from bokeh.io import show
from bokeh.models import CheckboxButtonGroup, CustomJS

LABELS = ["Option 1", "Option 2", "Option 3"]

checkbox_button_group = CheckboxButtonGroup(labels=LABELS, active=[0, 1])
checkbox_button_group.js_on_click(CustomJS(code="""
    console.log('checkbox_button_group: active=' + this.active, this.toString())
"""))

show(checkbox_button_group)
1
2
3
4
5
6
7
8
9
10
11
from bokeh.io import show
from bokeh.models import CheckboxGroup, CustomJS

LABELS = ["Option 1", "Option 2", "Option 3"]

checkbox_group = CheckboxGroup(labels=LABELS, active=[0, 1])
checkbox_group.js_on_click(CustomJS(code="""
    console.log('checkbox_group: active=' + this.active, this.toString())
"""))

show(checkbox_group)
1
2
temp = zip([1,2],['a','b'])
temp
1
<zip at 0x1cad5461880>
1
zip?
1
2
3
4
5
6
7
8
9
10
Init signature: zip(self, /, *args, **kwargs)
Docstring:   
zip(*iterables) --> zip object

Return a zip object whose .__next__() method returns a tuple where
the i-th element comes from the i-th iterable argument.  The .__next__()
method continues until the shortest iterable in the argument sequence
is exhausted and then it raises StopIteration.
Type:           type
Subclasses:   
1

1

1

1

1

1

1

1
2
3
4
5
import pandas as pd
import numpy as np

dates = pd.date_range('20130101', periods=6)
df = pd.DataFrame(np.random.randn(6, 4), index=dates, columns=list('ABCD'))
1
df
1
df[df > 0]
1
df[df['A'].isin(['a', 'b'])]
1

1

1
2
3
4
5
6
import matplotlib.pyplot as plt
import numpy as np
T = np.array([6, 7, 8, 9, 10, 11, 12])
power = np.array([1.53E+03, 5.92E+02, 2.04E+02, 7.24E+01, 2.72E+01, 1.10E+01, 4.70E+00])
plt.plot(T,power)
plt.show()

png

1
2
3
4
5
6
7
8
9
10
import matplotlib.pyplot as plt
import numpy as np
T = np.array([6, 7, 8, 9, 10, 11, 12])
power = np.array([1.53E+03, 5.92E+02, 2.04E+02, 7.24E+01, 2.72E+01, 1.10E+01, 4.70E+00])

from scipy.interpolate import spline
xnew = np.linspace(T.min(),T.max(),300) #300 represents number of points to make between T.min and T.max
power_smooth = spline(T,power,xnew)
plt.plot(xnew,power_smooth)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
---------------------------------------------------------------------------

ImportError                               Traceback (most recent call last)

<ipython-input-2-63c9a1a8c2a4> in <module>
      4 power = np.array([1.53E+03, 5.92E+02, 2.04E+02, 7.24E+01, 2.72E+01, 1.10E+01, 4.70E+00])
      5 
----> 6 from scipy.interpolate import spline
      7 xnew = np.linspace(T.min(),T.max(),300) #300 represents number of points to make between T.min and T.max
      8 power_smooth = spline(T,power,xnew)


ImportError: cannot import name 'spline' from 'scipy.interpolate' (E:\ProgramFiles\Anaconda\lib\site-packages\scipy\interpolate\__init__.py)
1

1

OS模块

  • os.getcwd()
  • os.name()
  • os.path()
  • os.listdir()
  • os.rename()
1
2
3
4
import os
print(os.getcwd())
print(os.name)
os.mkdir('./OS创建文件夹')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
E:\GitHubWP\JupyterWP
nt



---------------------------------------------------------------------------

FileExistsError                           Traceback (most recent call last)

<ipython-input-1-d2c10eafdc92> in <module>
      2 print(os.getcwd())
      3 print(os.name)
----> 4 os.mkdir('./OS创建文件夹')


FileExistsError: [WinError 183] 当文件已存在时，无法创建该文件。: './OS创建文件夹'
1
2
3
4
5
6
7
8
9
import os.path as op

# isdir(path)
# isfile(path)
# splitext(path) # 分开后缀
# join()
# exists() # 判断文件（夹）是否存在
# getsize() 
# getctime() # 获取创建时间,返回时间戳（1970-01-01，00:00:00）
1
2
results = op.abspath('./OS创建文件夹') # 返回绝对路径
results
1
'E:\\GitHubWP\\JupyterWP\\OS创建文件夹'
1
# dir(op)
1
2
result = op.getctime('./OS创建文件夹')
print(result)
1
1609938276.518
1
2
3
4
5
6
import time
print(result)
timearr = time.localtime(result)
print(timearr)
mytime = time.strftime('%Y-%m-%d %H:%M:%S',timearr)
print(mytime)
1
2
3
1609938276.518
time.struct_time(tm_year=2021, tm_mon=1, tm_mday=6, tm_hour=21, tm_min=4, tm_sec=36, tm_wday=2, tm_yday=6, tm_isdst=0)
2021-01-06 21:04:36
1
os.listdir(os.getcwd())[:5] # os.listdir会列出隐藏文件
1
2
3
4
5
['.ipynb_checkpoints',
 '0. python日积月累.ipynb',
 '2014-master',
 'alien_invasion',
 'All数据寻峰-面积-元素-含量表.xlsx']
1
2
3
4
5
# 批量修改文件名
# 1.获取指定的路径
# 2.获取路径下的文件名
# 3.修改文件名
# 4.批量？循环 or 函数！
1
2
3
4
5
6
7
8
9
10
import os
import os.path as op
def mkdir(dirname):
    path = op.join(os.getcwd(),dirname)
    folder = op.exists(path)
 
    if not folder:                   #判断是否存在文件夹如果不存在则创建为文件夹
        os.makedirs(path)            #makedirs 创建文件时如果路径不存在会创建这个路径

mkdir('testdir\子目录')
1

json模块

JSON（JavaScript Object Notation)是一种轻量级的数据交换格式，易于人阅读和编写，同也易于机器解析和生成，并有效地提升网络传输效率。

给字典加上引号，变成字符串，即json格式字符串。

1
2
3
4
5
6
7
8
9
10
11
12
json_str = '{"name":"张三"}'
json_str1 = '["name","张三"]'

# 将json格式的str转换成python的数据格式
import json
result = json.loads(json_str)
print(result,type(result))
result1 = json.loads(json_str1)
print(result1,type(result1))

result2 = json.dumps(result,ensure_ascii=False)
print(result2,type(result2))
1
2
3
{'name': '张三'} <class 'dict'>
['name', '张三'] <class 'list'>
{"name": "张三"} <class 'str'>
1

1

1

1

1

1

1

1

基础知识新知

函数新知

“代码有很多种坏味道，重复是最坏的一种”

  1. 默认参数
  2. 可变参数
  3. 同名函数：python没有函数重载的概念，后面的定义会覆盖前面的定义。
  • 用模块管理
  1. 模块中编写执行代码
    if __name__ == '__main__'
  2. 变量作用域

格式化输出的几种方式

1
2
a, b = 5, 10
print('%d * %d = %d' % (a, b, a * b))
1
5 * 10 = 50
1
2
a, b = 5, 10
print('{0} * {1} = {2}'.format(a, b, a * b))
1
5 * 10 = 50
1
2
a, b = 5, 10
print(f'{a} * {b} = {a * b}')
1
5 * 10 = 50

列表

1
[[1,2]]*5
1
[[1, 2], [1, 2], [1, 2], [1, 2], [1, 2]]

输入

1
a,b =map(int,input('输入a,b逗号隔开:').split(','))
1
输入a,b逗号隔开: 1,2
1
a
1
1
1
b
1
2

分形

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
import numpy as np
import matplotlib.pyplot as plt
def mandelbrot( h,w, maxit=20 ):
    """Returns an image of the Mandelbrot fractal of size (h,w)."""
    y,x = np.ogrid[ -1.4:1.4:h*1j, -2:0.8:w*1j ]
    c = x+y*1j
    z = c
    divtime = maxit + np.zeros(z.shape, dtype=int)

    for i in range(maxit):
        z = z**2 + c
        diverge = z*np.conj(z) > 2**2            # who is diverging
        div_now = diverge & (divtime==maxit)  # who is diverging now
        divtime[div_now] = i                  # note when
        z[diverge] = 2                        # avoid diverging too much

    return divtime
plt.imshow(mandelbrot(400,400))
#plt.savefig('scatter.eps',dpi=600,format='eps')
plt.show()

png

批量下载脚本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# 该例子来源于知乎

import urllib.request# url request
import re            # regular expression
import os            # dirs
import time

'''
url 下载网址
pattern 正则化的匹配关键词
Directory 下载目录
'''
def BatchDownload(url,pattern,Directory):
  
    # 拉动请求，模拟成浏览器去访问网站->跳过反爬虫机制
    headers = {'User-Agent', 'Mozilla/5.0 (Windows NT 6.1; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/64.0.3282.186 Safari/537.36'}
    opener = urllib.request.build_opener()
    opener.addheaders = [headers]
  
    # 获取网页内容
    content = opener.open(url).read().decode('utf8')
  
    # 构造正则表达式，从content中匹配关键词pattern
    raw_hrefs = re.findall(pattern, content, 0)
  
    # set函数消除重复元素
    hset = set(raw_hrefs)
       
    # 下载链接
    for href in hset:
        # 之所以if else 是为了区别只有一个链接的特别情况
        if(len(hset)>1):
            link = url + href[0]
            filename = os.path.join(Directory, href[0])
            print("正在下载",filename)
            urllib.request.urlretrieve(link, filename)
            print("成功下载！")
        else:
            link = url +href
            filename = os.path.join(Directory, href)
            print("正在下载",filename)
            urllib.request.urlretrieve(link, filename)
            print("成功下载！")
          
        # 无sleep间隔，网站认定这种行为是攻击，反反爬虫
        time.sleep(1)
      
#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/',
#              '(Storm-Data-Export-Format.docx)',
#              'E:\stormevents\csvfiles')
      
#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/',
#              '(Storm-Data-Export-Format.pdf)',
#              'E:\stormevents\csvfiles')
      
#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/',
#              '(StormEvents_details-ftp_v1.0_d(\d*)_c(\d*).csv.gz)',
#              'E:\stormevents\csvfiles')
      
#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/',
#              '(StormEvents_fatalities-ftp_v1.0_d(\d*)_c(\d*).csv.gz)',
#              'E:\stormevents\csvfiles')

#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/',
#              '(StormEvents_locations-ftp_v1.0_d(\d*)_c(\d*).csv.gz)',
#              'E:\stormevents\csvfiles')

#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/legacy/',
#              '(ugc_areas.csv)',
#              'E:\stormevents\csvfiles\legacy')

#BatchDownload('https://www1.ncdc.noaa.gov/pub/data/swdi/stormevents/csvfiles/',
#              '(ugc_areas.csv)',
#              'E:\stormevents\csvfiles')

3Blue1Brown 线性代数

  • 向量一般用列向量表示

  • 线性变换矩阵：由变换后的基向量组成

    • 逆时针旋转九十度：
    • 旋转矩阵：
    • 线性变换是操作空间的一种手段
    • 列线性相关，则降维
    • 将矩阵看做空间变换

  • 矩阵乘法与线性变换复合

    • 矩阵运算
    • 矩阵乘法不满足交换律
    • 矩阵乘法满足结合律（交换变换顺序）
  • 三维空间中的线性变换

  • 行列式

    • 行列式可以理解为变换前后的伸缩比例，复数表示改变了图形的定向（感觉上像是将空间翻转了：比如，j从i的左边变换到了右边）（三维空间中的定向：右手定则）
    • 只需要检验一个矩阵的行列式是否为零，就可以了解矩阵所代表的变换是否将空间压缩到更小的维度上
    • 计算公式：

bc项可理解为平行四边形法则中对角线的伸缩比例。

可用伸缩比例理解。

  • 逆矩阵、列空间、秩、零空间

    • 矩阵的用途：
      • 操作空间
      • 求解线性方程组
    • 线性方程组：求解一个向量经系数矩阵变换后得到已知的向量。
      • 系数矩阵的行列式为零
      • 系数矩阵的行列式不为零
        • 逆变换求解：逆矩阵
        • 理解矩阵与逆矩阵的乘法等于单位矩阵：恒等变换
        • 矩阵行列式不为零，存在逆矩阵。
        • 矩阵行列式为零，则矩阵变换将空间压缩到更低的维度上，没有逆矩阵（例，不能将一条线“解压缩”为一个平面）
        • 即便行列式为零，仍然可能存在解（例，已知向量正好在压缩后的直线上）
    • 秩
      • 秩代表着变换后空间的维数。
    • 列空间：矩阵的列所张成的空间
      • 秩是列空间的维数
      • 满秩
      • 零向量一定为被包含在列空间中（因为线性变换必须保持原点位置不变；线性变换是基于原点的旋转伸缩变换）
    • 零空间：
      • 核：变换后落在原点的向量的集合，称为矩阵的零空间或核。
      • 当已知向量为零向量时，零空间的向量就是对应的解。
    • 非方阵的情况
      • 2*矩阵的集合几何意义是将二维空间映射到三维空间上，该矩阵的列空间是三维空间中一个过原点的二维平面。
      • 其他变换：二维到一维，点积？

  • 点积与对偶性（daulity）

    • 几何解释：一个向量的长度乘以另一个向量的投影长度
    • 交换顺序，结果不变
    • 投影矩阵：与单位向量的点积可以理解为将向量投影到单位向量所在的直线上所得到的投影长度。
    • 点积：列向量转转置为行向量，即投影矩阵。

  • 叉积

    • 不严谨解释
      • 平行四边形的面积
      • 定向：i*j = 1;i在j的右侧为正，否则为负；夹角的方向定义为从前向量转到后向量的方向。
      • 二维向量可用行列式计算。（转置不改变行列式的值）
      • 向量垂直时，面积最大
    • 严谨解释
      • 真正的叉积是通过两个三维向量产生一个新的三维向量。
    • 混合积：
    • 对偶向量

  • 基变换：不同坐标系（不同视角）

    • 原点相同；基底不同
    • 对方的基底在我方坐标系下的坐标组成变换矩阵A.（该矩阵将我方坐标网格转换为对方网格：对同一向量而言，则相当于将对方语言翻译为我方语言）
    • 对方视角下的向量v,在我方看来则为Av.(类似将对方的语言翻译为我方的语言)
    • 该矩阵A取逆,记为invA，将对方网格转换为我方网格，对同一向量而言，相当于将我方语言转换为对方语言。invA v.

线性变换不同于坐标变换，线性变换是向量在同一坐标系下的伸缩转动；坐标变换是同一向量在不同坐标系下的坐标表示。

1
2
- 如何变换一个矩阵：对方基向量在我方表示为A.我方做了一个变换M,则其基向量跟着变换为MA，但这是我方的语言，在对方看来，应为 invAMA.
- invAMA暗示一种数学上的变换，中间的矩阵代表一种变换（例如，旋转90度），两边的矩阵代表视角上的变化。

  • 特征值与特征向量

    • 特征向量：在变换中留在其张成空间中的向量；变化仅有伸缩倍数，即特征值。特征值可以为负值。
    • 三维旋转：三维空间中的特征向量就是旋转轴；把一个三维旋转看成绕某个轴旋转一定角度，比考虑响应的3*3矩阵直观得多。（这种情况下特征值必须为1，因为只旋转，无伸缩）
    • 求解特征值与特征向量Mv=bv
      • 改变形式使右边为零：（M-bI）v=0
      • 右边为零向量，即被压缩，则行列式为0.可以求得特征值；特征值非实数表明没有特征向量
      • 再根据特征值，像求解线性方程组那样，可以求解特征向量。
    • 存在一个特征值，多个特征向量的情况，比如讲坐标拉伸2倍，则特征值为2，所有的向量都是特征向量。
    • 特殊情况：基向量是特征向量；对角矩阵；对角元是特征值
      • 对角矩阵乘法的计算
      • 用特征向量作为基，进行变换；换到特征向量作基底的视角下，该视角下的变换矩阵必然是对角的，且对角元为特征值。（前提，特征向量能张成全空间）；可应用于矩阵运算

  • 抽象向量空间

    • 行列式和特征向量不受所选坐标系的影响
      • 行列式表示缩放比例
      • 特征向量表示在变换中留在它所张成的空间中的向量
      • 你可以自由选取坐标系，这并不会改变它们最根本的值。
    • 空间性
    • 函数式另一种向量？
      • 向量的所有特征均可应用于函数
        • 相加
        • 数乘
        • 线性变换：符合函数、求导
      • 线性算子：一个函数式线性的？
        • 线性的严格定义：可加性；成比例（一阶齐次）
        • 线性变换保持向量加法运算和数乘运算
        • 求导是线性运算
        • 基函数
        • 求导算符的矩阵：选择基函数，建立矩阵
    • 向量是什么？向量的形式不重要，只要满足向量定义的一组格则，它可以是任何东西。

  • 克莱姆法则，几何解释

    • 计算线性方程组

      • 高斯消元法
      • 克莱姆法则

      不改变点积的变换称为正交变换，例如旋转矩阵

pyecharts

1
2
from pyecharts.globals import CurrentConfig, NotebookType
CurrentConfig.NOTEBOOK_TYPE = NotebookType.JUPYTER_LAB
1
2
3
4
5
6
7
8
from pyecharts.charts import Bar

bar = Bar()
bar.add_xaxis(["衬衫", "羊毛衫", "雪纺衫", "裤子", "高跟鞋", "袜子"])
bar.add_yaxis("商家A", [5, 20, 36, 10, 75, 90])
# render 会生成本地 HTML 文件，默认会在当前目录生成 render.html 文件
# 也可以传入路径参数，如 bar.render("mycharts.html")
bar.render_notebook()

<!DOCTYPE html>

Geoprocessing with python

1
from osgeo import ogr
1
2
driver = ogr.GetDriverByName('geojson')
print(driver)
1
<osgeo.ogr.Driver; proxy of <Swig Object of type 'OGRDriverShadow *' at 0x00000228D2B5C330> >
1
2
import ospybook as pb
pb.print.drivers()
1
2
3
4
5
6
7
8
9
10
---------------------------------------------------------------------------

ModuleNotFoundError                       Traceback (most recent call last)

<ipython-input-3-4937f463c8c4> in <module>
----> 1 import ospybook as pb
      2 pb.print.drivers()


ModuleNotFoundError: No module named 'ospybook'
1
# %pip install ospybook
1
2
3
4
5
Note: you may need to restart the kernel to use updated packages.


ERROR: Could not find a version that satisfies the requirement ospybook (from versions: none)
ERROR: No matching distribution found for ospybook
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
"""
Definition of colour schemes for lines and maps that also work for colour-blind
people. See https://personal.sron.nl/~pault/ for background information and
best usage of the schemes.

Copyright (c) 2019, Paul Tol
All rights reserved.

License:  Standard 3-clause BSD
"""
import numpy as np
from matplotlib.colors import LinearSegmentedColormap, to_rgba_array


def discretemap(colormap, hexclrs):
    """
    Produce a colormap from a list of discrete colors without interpolation.
    """
    clrs = to_rgba_array(hexclrs)
    clrs = np.vstack([clrs[0], clrs, clrs[-1]])
    cdict = {}
    for ki, key in enumerate(('red','green','blue')):
        cdict[key] = [ (i/(len(clrs)-2.), clrs[i, ki], clrs[i+1, ki]) for i in range(len(clrs)-1) ]
    return LinearSegmentedColormap(colormap, cdict)


class TOLcmaps(object):
    """
    Class TOLcmaps definition.
    """
    def __init__(self):
        """
        """
        self.cmap = None
        self.cname = None
        self.namelist = (
            'sunset_discrete', 'sunset', 'BuRd_discrete', 'BuRd',
            'PRGn_discrete', 'PRGn', 'YlOrBr_discrete', 'YlOrBr', 'WhOrBr',
            'iridescent', 'rainbow_PuRd', 'rainbow_PuBr', 'rainbow_WhRd',
            'rainbow_WhBr', 'rainbow_discrete')

        self.funcdict = dict(
            zip(self.namelist,
                (self.__sunset_discrete, self.__sunset, self.__BuRd_discrete,
                self.__BuRd, self.__PRGn_discrete, self.__PRGn,
                self.__YlOrBr_discrete, self.__YlOrBr, self.__WhOrBr,
                self.__iridescent, self.__rainbow_PuRd, self.__rainbow_PuBr,
                self.__rainbow_WhRd, self.__rainbow_WhBr,
                self.__rainbow_discrete)))

    def __sunset_discrete(self):
        """
        Define colormap 'sunset_discrete'.
        """
        clrs = ['#364B9A', '#4A7BB7', '#6EA6CD', '#98CAE1', '#C2E4EF',
                '#EAECCC', '#FEDA8B', '#FDB366', '#F67E4B', '#DD3D2D',
                '#A50026']
        self.cmap = discretemap(self.cname, clrs)
        self.cmap.set_bad('#FFFFFF')

    def __sunset(self):
        """
        Define colormap 'sunset'.
        """
        clrs = ['#364B9A', '#4A7BB7', '#6EA6CD', '#98CAE1', '#C2E4EF',
                '#EAECCC', '#FEDA8B', '#FDB366', '#F67E4B', '#DD3D2D',
                '#A50026']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#FFFFFF')

    def __BuRd_discrete(self):
        """
        Define colormap 'BuRd_discrete'.
        """
        clrs = ['#2166AC', '#4393C3', '#92C5DE', '#D1E5F0', '#F7F7F7',
                '#FDDBC7', '#F4A582', '#D6604D', '#B2182B']
        self.cmap = discretemap(self.cname, clrs)
        self.cmap.set_bad('#FFEE99')

    def __BuRd(self):
        """
        Define colormap 'BuRd'.
        """
        clrs = ['#2166AC', '#4393C3', '#92C5DE', '#D1E5F0', '#F7F7F7',
                '#FDDBC7', '#F4A582', '#D6604D', '#B2182B']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#FFEE99')

    def __PRGn_discrete(self):
        """
        Define colormap 'PRGn_discrete'.
        """
        clrs = ['#762A83', '#9970AB', '#C2A5CF', '#E7D4E8', '#F7F7F7',
                '#D9F0D3', '#ACD39E', '#5AAE61', '#1B7837']
        self.cmap = discretemap(self.cname, clrs)
        self.cmap.set_bad('#FFEE99')

    def __PRGn(self):
        """
        Define colormap 'PRGn'.
        """
        clrs = ['#762A83', '#9970AB', '#C2A5CF', '#E7D4E8', '#F7F7F7',
                '#D9F0D3', '#ACD39E', '#5AAE61', '#1B7837']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#FFEE99')

    def __YlOrBr_discrete(self):
        """
        Define colormap 'YlOrBr_discrete'.
        """
        clrs = ['#FFFFE5', '#FFF7BC', '#FEE391', '#FEC44F', '#FB9A29',
                '#EC7014', '#CC4C02', '#993404', '#662506']
        self.cmap = discretemap(self.cname, clrs)
        self.cmap.set_bad('#888888')

    def __YlOrBr(self):
        """
        Define colormap 'YlOrBr'.
        """
        clrs = ['#FFFFE5', '#FFF7BC', '#FEE391', '#FEC44F', '#FB9A29',
                '#EC7014', '#CC4C02', '#993404', '#662506']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#888888')

    def __WhOrBr(self):
        """
        Define colormap 'WhOrBr'.
        """
        clrs = ['#FFFFFF', '#FFF7BC', '#FEE391', '#FEC44F', '#FB9A29',
                '#EC7014', '#CC4C02', '#993404', '#662506']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#888888')

    def __iridescent(self):
        """
        Define colormap 'iridescent'.
        """
        clrs = ['#FEFBE9', '#FCF7D5', '#F5F3C1', '#EAF0B5', '#DDECBF',
                '#D0E7CA', '#C2E3D2', '#B5DDD8', '#A8D8DC', '#9BD2E1',
                '#8DCBE4', '#81C4E7', '#7BBCE7', '#7EB2E4', '#88A5DD',
                '#9398D2', '#9B8AC4', '#9D7DB2', '#9A709E', '#906388',
                '#805770', '#684957', '#46353A']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#999999')

    def __rainbow_PuRd(self):
        """
        Define colormap 'rainbow_PuRd'.
        """
        clrs = ['#6F4C9B', '#6059A9', '#5568B8', '#4E79C5', '#4D8AC6',
                '#4E96BC', '#549EB3', '#59A5A9', '#60AB9E', '#69B190',
                '#77B77D', '#8CBC68', '#A6BE54', '#BEBC48', '#D1B541',
                '#DDAA3C', '#E49C39', '#E78C35', '#E67932', '#E4632D',
                '#DF4828', '#DA2222']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#FFFFFF')

    def __rainbow_PuBr(self):
        """
        Define colormap 'rainbow_PuBr'.
        """
        clrs = ['#6F4C9B', '#6059A9', '#5568B8', '#4E79C5', '#4D8AC6',
                '#4E96BC', '#549EB3', '#59A5A9', '#60AB9E', '#69B190',
                '#77B77D', '#8CBC68', '#A6BE54', '#BEBC48', '#D1B541',
                '#DDAA3C', '#E49C39', '#E78C35', '#E67932', '#E4632D',
                '#DF4828', '#DA2222', '#B8221E', '#95211B', '#721E17',
                '#521A13']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#FFFFFF')

    def __rainbow_WhRd(self):
        """
        Define colormap 'rainbow_WhRd'.
        """
        clrs = ['#E8ECFB', '#DDD8EF', '#D1C1E1', '#C3A8D1', '#B58FC2',
                '#A778B4', '#9B62A7', '#8C4E99', '#6F4C9B', '#6059A9',
                '#5568B8', '#4E79C5', '#4D8AC6', '#4E96BC', '#549EB3',
                '#59A5A9', '#60AB9E', '#69B190', '#77B77D', '#8CBC68',
                '#A6BE54', '#BEBC48', '#D1B541', '#DDAA3C', '#E49C39',
                '#E78C35', '#E67932', '#E4632D', '#DF4828', '#DA2222']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#666666')

    def __rainbow_WhBr(self):
        """
        Define colormap 'rainbow_WhBr'.
        """
        clrs = ['#E8ECFB', '#DDD8EF', '#D1C1E1', '#C3A8D1', '#B58FC2',
                '#A778B4', '#9B62A7', '#8C4E99', '#6F4C9B', '#6059A9',
                '#5568B8', '#4E79C5', '#4D8AC6', '#4E96BC', '#549EB3',
                '#59A5A9', '#60AB9E', '#69B190', '#77B77D', '#8CBC68',
                '#A6BE54', '#BEBC48', '#D1B541', '#DDAA3C', '#E49C39',
                '#E78C35', '#E67932', '#E4632D', '#DF4828', '#DA2222',
                '#B8221E', '#95211B', '#721E17', '#521A13']
        self.cmap = LinearSegmentedColormap.from_list(self.cname, clrs)
        self.cmap.set_bad('#666666')

    def __rainbow_discrete(self, lut=None):
        """
        Define colormap 'rainbow_discrete'.
        """
        clrs = ['#E8ECFB', '#D9CCE3', '#D1BBD7', '#CAACCB', '#BA8DB4',
                '#AE76A3', '#AA6F9E', '#994F88', '#882E72', '#1965B0',
                '#437DBF', '#5289C7', '#6195CF', '#7BAFDE', '#4EB265',
                '#90C987', '#CAE0AB', '#F7F056', '#F7CB45', '#F6C141',
                '#F4A736', '#F1932D', '#EE8026', '#E8601C', '#E65518',
                '#DC050C', '#A5170E', '#72190E', '#42150A']
        indexes = [[9], [9, 25], [9, 17, 25], [9, 14, 17, 25], [9, 13, 14, 17,
                25], [9, 13, 14, 16, 17, 25], [8, 9, 13, 14, 16, 17, 25], [8,
                9, 13, 14, 16, 17, 22, 25], [8, 9, 13, 14, 16, 17, 22, 25, 27],
                [8, 9, 13, 14, 16, 17, 20, 23, 25, 27], [8, 9, 11, 13, 14, 16,
                17, 20, 23, 25, 27], [2, 5, 8, 9, 11, 13, 14, 16, 17, 20, 23,
                25], [2, 5, 8, 9, 11, 13, 14, 15, 16, 17, 20, 23, 25], [2, 5,
                8, 9, 11, 13, 14, 15, 16, 17, 19, 21, 23, 25], [2, 5, 8, 9, 11,
                13, 14, 15, 16, 17, 19, 21, 23, 25, 27], [2, 4, 6, 8, 9, 11,
                13, 14, 15, 16, 17, 19, 21, 23, 25, 27], [2, 4, 6, 7, 8, 9, 11,
                13, 14, 15, 16, 17, 19, 21, 23, 25, 27], [2, 4, 6, 7, 8, 9, 11,
                13, 14, 15, 16, 17, 19, 21, 23, 25, 26, 27], [1, 3, 4, 6, 7, 8,
                9, 11, 13, 14, 15, 16, 17, 19, 21, 23, 25, 26, 27], [1, 3, 4,
                6, 7, 8, 9, 10, 12, 13, 14, 15, 16, 17, 19, 21, 23, 25, 26,
                27], [1, 3, 4, 6, 7, 8, 9, 10, 12, 13, 14, 15, 16, 17, 18, 20,
                22, 24, 25, 26, 27], [1, 3, 4, 6, 7, 8, 9, 10, 12, 13, 14, 15,
                16, 17, 18, 20, 22, 24, 25, 26, 27, 28], [0, 1, 3, 4, 6, 7, 8,
                9, 10, 12, 13, 14, 15, 16, 17, 18, 20, 22, 24, 25, 26, 27, 28]]
        if lut == None or lut < 1 or lut > 23:
            lut = 22
        self.cmap = discretemap(self.cname, [ clrs[i] for i in indexes[lut-1] ])
        if lut == 23:
            self.cmap.set_bad('#777777')
        else:
            self.cmap.set_bad('#FFFFFF')

    def show(self):
        """
        List names of defined colormaps.
        """
        print(' '.join(repr(n) for n in self.namelist))

    def get(self, cname='rainbow_PuRd', lut=None):
        """
        Return requested colormap, default is 'rainbow_PuRd'.
        """
        self.cname = cname
        if cname == 'rainbow_discrete':
            self.__rainbow_discrete(lut)
        else:
            self.funcdict[cname]()
        return self.cmap


def tol_cmap(colormap=None, lut=None):
    """
    Continuous and discrete color sets for ordered data.
  
    Return a matplotlib colormap.
    Parameter lut is ignored for all colormaps except 'rainbow_discrete'.
    """
    obj = TOLcmaps()
    if colormap == None:
        return obj.namelist
    if colormap not in obj.namelist:
        colormap = 'rainbow_PuRd'
        print('*** Warning: requested colormap not defined,',
              'known colormaps are {}.'.format(obj.namelist),
              'Using {}.'.format(colormap))
    return obj.get(colormap, lut)


def tol_cset(colorset=None):
    """
    Discrete color sets for qualitative data.

    Define a namedtuple instance with the colors.
    Examples for: cset = tol_cset(<scheme>)
      - cset.red and cset[1] give the same color (in default 'bright' colorset)
      - cset._fields gives a tuple with all color names
      - list(cset) gives a list with all colors
    """
    from collections import namedtuple
  
    namelist = ('bright', 'high-contrast', 'vibrant', 'muted', 'light')
    if colorset == None:
        return namelist
    if colorset not in namelist:
        colorset = 'bright'
        print('*** Warning: requested colorset not defined,',
              'known colorsets are {}.'.format(namelist),
              'Using {}.'.format(colorset)) 

    if colorset == 'bright':
        cset = namedtuple('Bcset',
                    'blue red green yellow cyan purple grey black')
        return cset('#4477AA', '#EE6677', '#228833', '#CCBB44', '#66CCEE',
                    '#AA3377', '#BBBBBB', '#000000')
  
    if colorset == 'high-contrast':
        cset = namedtuple('Hcset',
                    'blue yellow red black')
        return cset('#004488', '#DDAA33', '#BB5566', '#000000')

    if colorset == 'vibrant':
        cset = namedtuple('Vcset',
                    'orange blue cyan magenta red teal grey black')
        return cset('#EE7733', '#0077BB', '#33BBEE', '#EE3377', '#CC3311',
                    '#009988', '#BBBBBB', '#000000')

    if colorset == 'muted':
        cset = namedtuple('Mcset',
                    'rose indigo sand green cyan wine teal olive purple pale_grey black')
        return cset('#CC6677', '#332288', '#DDCC77', '#117733', '#88CCEE',
                    '#882255', '#44AA99', '#999933', '#AA4499', '#DDDDDD',
                    '#000000')

    if colorset == 'light':
        cset = namedtuple('Lcset',
                    'light_blue orange light_yellow pink light_cyan mint pear olive pale_grey black')
        return cset('#77AADD', '#EE8866', '#EEDD88', '#FFAABB', '#99DDFF',
                    '#44BB99', '#BBCC33', '#AAAA00', '#DDDDDD', '#000000')


def main():

    from matplotlib import pyplot as plt

    # Change default colorset (for lines) and colormap (for maps).
#    plt.rc('axes', prop_cycle=plt.cycler('color', list(tol_cset('bright'))))
#    plt.cm.register_cmap('rainbow_PuRd', tol_cmap('rainbow_PuRd'))
#    plt.rc('image', cmap='rainbow_PuRd')

    # Show colorsets tol_cset(<scheme>).
    schemes = tol_cset()
    fig, axes = plt.subplots(ncols=len(schemes))
    fig.subplots_adjust(top=0.92, bottom=0.02, left=0.02, right=0.92)
    for ax, scheme in zip(axes, schemes):
        cset = tol_cset(scheme)
        names = cset._fields
        colors = list(cset)
        for name, color in zip(names, colors):
            ax.scatter([], [], c=color, s=80, label=name)
        ax.set_axis_off()
        ax.legend(loc=2)
        ax.set_title(scheme)
    plt.show()

    # Show colormaps tol_cmap(<scheme>). 
    schemes = tol_cmap()
    gradient = np.linspace(0, 1, 256)
    gradient = np.vstack((gradient, gradient))
    fig, axes = plt.subplots(nrows=len(schemes))
    fig.subplots_adjust(top=0.98, bottom=0.02, left=0.2, right=0.99)
    for ax, scheme in zip(axes, schemes):
        pos = list(ax.get_position().bounds)
        ax.set_axis_off()
        ax.imshow(gradient, aspect=4, cmap=tol_cmap(scheme))
        fig.text(pos[0] - 0.01, pos[1] + pos[3]/2., scheme, va='center', ha='right', fontsize=10)
    plt.show()

    # Show colormaps tol_cmap('rainbow_discrete', <lut>). 
    gradient = np.linspace(0, 1, 256)
    gradient = np.vstack((gradient, gradient))
    fig, axes = plt.subplots(nrows=23)
    fig.subplots_adjust(top=0.98, bottom=0.02, left=0.25, right=0.99)
    for lut, ax in enumerate(axes, start=1):
        pos = list(ax.get_position().bounds)
        ax.set_axis_off()
        ax.imshow(gradient, aspect=4, cmap=tol_cmap('rainbow_discrete', lut))
        fig.text(pos[0] - 0.01, pos[1] + pos[3]/2., 'rainbow_discrete, ' + str(lut), va='center', ha='right', fontsize=10)
    plt.show()


if __name__ == '__main__':
    main()

png

png

png

正态分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

mu,sigma = 0,0.1
s = np.random.normal(mu,sigma,1000)

mean = abs(mu - np.mean(s))
print('mean:',mean)

var = abs(sigma - np.std(s,ddof=1))
print('var:',var)

count,bins,ignored = plt.hist(s,30,density=True)
plt.plot(bins,1/(sigma*np.sqrt(2*np.pi))*np.exp(-(bins-mu)**2/(2*sigma**2)),lw=2,c='r')
plt.show()
1
2
mean: 0.005425582125086683
var: 0.0023657053391205884

png

1

差分法求解微分方程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np
import matplotlib.pyplot as plt
rc = 2.0 #设置常数
dt = 0.5 #设置步长
n = 1000 #设置分割段数
t = 0.0 #设置初始时间
q = 1.0 #设置初始电量

#先定义三个空列表
qt=[] #用来盛放差分得到的q值
qt0=[] #用来盛放解析得到的q值
time = [] #用来盛放时间值

for i in range(n):
    t = t + dt
    q1 = q - q*dt/rc #qn+1的近似值
    q = q - 0.5*(q1*dt/rc + q*dt/rc) #差分递推关系
    q0 = np.exp(-t/rc) #解析关系
    qt.append(q) #差分得到的q值列表
    qt0.append(q0) #解析得到的q值列表
    time.append(t) #时间列表

plt.plot(time,qt,'o',label='Euler-Modify') #差分得到的电量随时间的变化
plt.plot(time,qt0,'r-',label='Analytical') #解析得到的电量随时间的变化
plt.xlabel('time')
plt.ylabel('charge')
plt.xlim(0,20)
plt.ylim(-0.2,1.0)
plt.legend(loc='upper right')
plt.show()

png

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np
import matplotlib.pyplot as plt

h = 0.1#空间步长
N =30#空间步数
dt = 0.0001#时间步长
M = 10000#时间的步数
A = dt/(h**2) #lambda*tau/h^2
U = np.zeros([N+1,M+1])#建立二维空数组
Space = np.arange(0,(N+1)*h,h)#建立空间等差数列，从0到3，公差是h

#边界条件
for k in np.arange(0,M+1):
    U[0,k] = 0.0
    U[N,k] = 0.0

#初始条件
for i in np.arange(0,N):
    U[i,0]=4*i*h*(3-i*h)

#递推关系
for k in np.arange(0,M):
    for i in np.arange(1,N):
        U[i,k+1]=A*U[i+1,k]+(1-2*A)*U[i,k]+A*U[i-1,k]
1
2
3
4
5
6
7
8
9
10
11
12
#不同时刻的温度随空间坐标的变化
plt.plot(Space,U[:,0], 'g-', label='t=0',linewidth=1.0)
plt.plot(Space,U[:,3000], 'b-', label='t=3/10',linewidth=1.0)
plt.plot(Space,U[:,6000], 'k-', label='t=6/10',linewidth=1.0)
plt.plot(Space,U[:,9000], 'r-', label='t=9/10',linewidth=1.0)
plt.plot(Space,U[:,10000], 'y-', label='t=1',linewidth=1.0)
plt.ylabel('u(x,t)', fontsize=20)
plt.xlabel('x', fontsize=20)
plt.xlim(0,3)
plt.ylim(-2,10)
plt.legend(loc='upper right')
plt.show()

png

1
2
3
4
5
6
7
#温度等高线随时空坐标的变化，温度越高，颜色越偏红
extent = [0,1,0,3]#时间和空间的取值范围
levels = np.arange(0,10,0.1)#温度等高线的变化范围0-10，变化间隔为0.1
plt.contourf(U,levels,origin='lower',extent=extent,cmap=plt.cm.jet)
plt.ylabel('x', fontsize=20)
plt.xlabel('t', fontsize=20)
plt.show()

png

最小二乘

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
from scipy.optimize import least_squares

def model(x, u):
    return x[0] * (u ** 2 + x[1] * u) / (u ** 2 + x[2] * u + x[3])

def fun(x, u, y):
    return model(x, u) - y

def jac(x, u, y):
    J = np.empty((u.size, x.size))
    den = u ** 2 + x[2] * u + x[3]
    num = u ** 2 + x[1] * u
    J[:, 0] = num / den
    J[:, 1] = x[0] * u / den
    J[:, 2] = -x[0] * num * u / den ** 2
    J[:, 3] = -x[0] * num / den ** 2
    return J

u = np.array([4.0, 2.0, 1.0, 5.0e-1, 2.5e-1, 1.67e-1, 1.25e-1, 1.0e-1,
              8.33e-2, 7.14e-2, 6.25e-2])
y = np.array([1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2,
              4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2])
x0 = np.array([2.5, 3.9, 4.15, 3.9])
res = least_squares(fun, x0, jac=jac, bounds=(0, 100), args=(u, y), verbose=2)



res.x


import matplotlib.pyplot as plt
u_test = np.linspace(0, 5)
y_test = model(res.x, u_test)
plt.plot(u, y, 'o', markersize=4, label='data')
plt.plot(u_test, y_test, label='fitted model')
plt.xlabel("u")
plt.ylabel("y")
plt.legend(loc='lower right')
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   
       0              1         4.4383e+00                                    1.04e+02  
       1              2         9.6619e-01      3.47e+00       5.94e+00       1.37e+01  
       2              3         1.8466e-01      7.82e-01       8.49e+00       1.87e+00  
       3              4         2.9667e-02      1.55e-01       9.68e+00       2.70e-01  
       4              5         5.8265e-03      2.38e-02       8.95e+00       3.72e-02  
       5              6         3.4795e-03      2.35e-03       5.49e+00       3.57e-02  
       6              9         3.0037e-03      4.76e-04       8.95e-01       4.42e-02  
       7             10         2.4853e-03      5.18e-04       2.18e+00       5.05e-01  
       8             12         2.1802e-03      3.05e-04       7.84e-01       9.54e-02  
       9             13         1.9495e-03      2.31e-04       1.56e+00       3.77e-01  
      10             14         1.8691e-03      8.04e-05       2.60e+00       1.14e+00  
      11             15         1.5814e-03      2.88e-04       6.36e-01       8.71e-03  
      12             16         1.4525e-03      1.29e-04       1.46e+00       6.88e-02  
      13             17         1.3709e-03      8.16e-05       3.02e+00       9.07e-01  
      14             18         1.2146e-03      1.56e-04       2.84e+00       6.88e-01  
      15             20         1.1383e-03      7.62e-05       1.40e+00       1.21e-01  
      16             21         1.0924e-03      4.59e-05       2.89e+00       6.50e-01  
      17             22         1.0322e-03      6.02e-05       2.80e+00       5.30e-01  
      18             24         1.0004e-03      3.18e-05       1.37e+00       1.01e-01  
      19             25         9.7276e-04      2.76e-05       2.71e+00       5.07e-01  
      20             26         9.4014e-04      3.26e-05       2.63e+00       4.04e-01  
      21             28         9.2354e-04      1.66e-05       1.28e+00       8.28e-02  
      22             29         9.0343e-04      2.01e-05       2.53e+00       3.95e-01  
      23             30         8.9155e-04      1.19e-05       4.82e+00       1.48e+00  
      24             31         8.3675e-04      5.48e-05       4.37e+00       9.37e-01  
      25             33         8.1455e-04      2.22e-05       1.76e+00       1.22e-01  
      26             34         7.9193e-04      2.26e-05       3.41e+00       7.37e-01  
      27             35         7.8059e-04      1.13e-05       6.50e+00       3.22e+00  
      28             36         7.1291e-04      6.77e-05       2.43e+00       1.00e-01  
      29             37         6.6760e-04      4.53e-05       5.32e+00       2.52e+00  
      30             38         6.1812e-04      4.95e-05       2.96e+00       6.26e-01  
      31             40         5.9431e-04      2.38e-05       1.52e+00       1.42e-01  
      32             41         5.4701e-04      4.73e-05       2.92e+00       8.77e-01  
      33             42         5.0317e-04      4.38e-05       4.21e+00       4.81e+00  
      34             43         4.7950e-04      2.37e-05       7.87e-01       6.81e-05  
      35             44         4.7049e-04      9.01e-06       8.77e+00       1.23e-02  
      36             45         4.6289e-04      7.60e-06       2.97e-01       1.55e-04  
      37             47         4.6277e-04      1.20e-07       2.88e+00       9.10e-01  
      38             48         4.6212e-04      6.53e-07       1.44e+00       2.69e-01  
      39             49         4.6196e-04      1.55e-07       2.82e+00       9.28e-01  
      40             50         4.6117e-04      7.87e-07       1.52e+00       3.13e-01  
      41             51         4.6114e-04      3.49e-08       2.93e+00       1.09e+00  
      42             52         4.6015e-04      9.92e-07       7.46e-01       8.75e-02  
      43             53         4.5976e-04      3.89e-07       1.50e+00       3.05e-01  
      44             55         4.5946e-04      3.01e-07       7.61e-01       8.79e-02  
      45             56         4.5903e-04      4.22e-07       1.53e+00       3.47e-01  
      46             58         4.5867e-04      3.61e-07       7.73e-01       9.77e-02  
      47             59         4.5820e-04      4.68e-07       1.56e+00       3.80e-01  
      48             61         4.5777e-04      4.30e-07       7.85e-01       1.07e-01  
      49             62         4.5725e-04      5.21e-07       1.58e+00       4.16e-01  
      50             64         4.5674e-04      5.18e-07       7.95e-01       1.18e-01  
      51             65         4.5615e-04      5.81e-07       1.60e+00       4.56e-01  
      52             67         4.5552e-04      6.32e-07       8.04e-01       1.31e-01  
      53             68         4.5488e-04      6.47e-07       1.61e+00       5.01e-01  
      54             70         4.5409e-04      7.82e-07       8.10e-01       1.44e-01  
      55             71         4.5337e-04      7.19e-07       1.62e+00       5.50e-01  
      56             72         4.5231e-04      1.06e-06       1.64e+00       6.29e-01  
      57             74         4.5113e-04      1.18e-06       7.19e-01       1.29e-01  
      58             75         4.5019e-04      9.46e-07       1.44e+00       4.98e-01  
      59             77         4.4908e-04      1.11e-06       7.22e-01       1.45e-01  
      60             78         4.4803e-04      1.05e-06       1.44e+00       5.53e-01  
      61             79         4.4649e-04      1.54e-06       1.44e+00       6.29e-01  
      62             81         4.4476e-04      1.73e-06       6.52e-01       1.37e-01  
      63             82         4.4334e-04      1.41e-06       1.30e+00       5.29e-01  
      64             84         4.4161e-04      1.74e-06       6.40e-01       1.49e-01  
      65             85         4.3999e-04      1.62e-06       1.27e+00       5.63e-01  
      66             86         4.3754e-04      2.45e-06       1.24e+00       6.16e-01  
      67             88         4.3507e-04      2.48e-06       5.08e-01       1.09e-01  
      68             89         4.3277e-04      2.29e-06       1.02e+00       4.35e-01  
      69             91         4.3059e-04      2.18e-06       4.84e-01       1.11e-01  
      70             92         4.2791e-04      2.68e-06       9.60e-01       4.35e-01  
      71             94         4.2539e-04      2.52e-06       4.42e-01       1.02e-01  
      72             95         4.2217e-04      3.22e-06       8.80e-01       4.07e-01  
      73             97         4.1937e-04      2.80e-06       3.97e-01       8.97e-02  
      74             98         4.1553e-04      3.84e-06       7.97e-01       3.74e-01  
      75             100        4.1245e-04      3.07e-06       3.53e-01       7.69e-02  
      76             101        4.0791e-04      4.55e-06       7.17e-01       3.39e-01  
      77             102        4.0609e-04      1.81e-06       1.33e+00       1.28e+00  
      78             103        3.9204e-04      1.41e-05       2.07e-01       6.47e-03  
      79             104        3.8550e-04      6.54e-06       5.50e-01       2.61e-01  
      80             105        3.7640e-04      9.10e-06       1.01e+00       9.69e-01  
      81             106        3.6016e-04      1.62e-05       3.44e-01       7.06e-02  
      82             108        3.5059e-04      9.57e-06       4.10e-01       2.10e-01  
      83             109        3.3248e-04      1.81e-05       7.43e-01       7.37e-01  
      84             110        3.1288e-04      1.96e-05       3.00e-01       7.35e-02  
      85             112        2.9892e-04      1.40e-05       2.97e-01       1.72e-01  
      86             113        2.6939e-04      2.95e-05       5.32e-01       5.82e-01  
      87             114        2.4550e-04      2.39e-05       2.12e-01       4.76e-02  
      88             116        2.2954e-04      1.60e-05       1.75e-01       1.09e-01  
      89             117        1.9634e-04      3.32e-05       3.11e-01       3.45e-01  
      90             118        1.7589e-04      2.04e-05       1.36e-01       3.46e-02  
      91             119        1.5868e-04      1.72e-05       2.53e-01       4.30e-01  
      92             120        1.5380e-04      4.88e-06       9.94e-03       2.48e-03  
      93             121        1.5376e-04      4.34e-08       6.30e-03       9.00e-04  
      94             122        1.5375e-04      3.43e-09       4.28e-03       4.76e-04  
      95             123        1.5375e-04      1.17e-09       2.57e-03       8.19e-05  
      96             124        1.5375e-04      4.50e-10       1.60e-03       2.06e-04  
      97             125        1.5375e-04      1.74e-10       1.00e-03       1.22e-05  
      98             126        1.5375e-04      6.80e-11       6.26e-04       8.48e-05  
      99             127        1.5375e-04      2.66e-11       3.92e-04       1.88e-06  
      100            128        1.5375e-04      1.05e-11       2.46e-04       3.40e-05  
      101            129        1.5375e-04      4.12e-12       1.54e-04       2.91e-07  
      102            130        1.5375e-04      1.62e-12       9.70e-05       1.35e-05  
      103            131        1.5375e-04      6.40e-13       6.09e-05       4.52e-08  
`ftol` termination condition is satisfied.
Function evaluations 131, initial cost 4.4383e+00, final cost 1.5375e-04, first-order optimality 4.52e-08.

png

1
res
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
active_mask: array([0, 0, 0, 0])
       cost: 0.0001537528023415085
        fun: array([-1.30664586e-03, -1.87579954e-03,  8.91958340e-03, -1.11095297e-02,
       8.35200541e-03, -1.73677561e-04,  2.58387522e-05,  1.26262347e-03,
      -3.52355168e-03,  6.16759030e-04, -3.88872553e-03])
       grad: array([-4.52564650e-10, -1.44529678e-10,  1.55211428e-09,  8.54196275e-08])
        jac: array([[ 1.00823291,  0.04638017, -0.04676201, -0.0116905 ],
      [ 1.0000944 ,  0.08799521, -0.08800352, -0.04400176],
      [ 0.94613022,  0.15312606, -0.14487719, -0.14487719],
      [ 0.77222944,  0.21537647, -0.16632005, -0.33264011],
      [ 0.4810639 ,  0.21017745, -0.10110878, -0.40443513],
      [ 0.32429656,  0.17450668, -0.05659192, -0.33887375],
      [ 0.23664118,  0.14424711, -0.03413481, -0.27307844],
      [ 0.18392905,  0.12173779, -0.02239112, -0.22391116],
      [ 0.14925079,  0.10479279, -0.01564041, -0.18775999],
      [ 0.12508303,  0.09180226, -0.01148291, -0.160825  ],
      [ 0.10742028,  0.08160364, -0.00876589, -0.14025417]])
    message: '`ftol` termination condition is satisfied.'
       nfev: 131
       njev: 104
 optimality: 4.516920778943166e-08
     status: 2
    success: True
          x: array([0.192806  , 0.19130332, 0.12306046, 0.13607205])
1
res.x
1
array([0.192806  , 0.19130332, 0.12306046, 0.13607205])
1
2
3
4
5
from scipy.optimize import leastsq
def func(x):
    return 2*(x-5)**2+x**3 -5*x -50
a = leastsq(func, 20)
print(a,func(a[0]))
1
(array([4.09901951]), 2) [-1.42108547e-14]
1
2
b = leastsq(func,-20,full_output=True)
b
1
2
3
4
5
6
7
8
9
(array([-6.09901951]),
 array([[0.00025849]]),
 {'fvec': array([3.55271368e-14]),
  'nfev': 17,
  'fjac': array([[-62.19803813]]),
  'ipvt': array([1], dtype=int32),
  'qtf': array([1.47077363e-06])},
 'The relative error between two consecutive iterates is at most 0.000000',
 2)
1
func(a[0])
1
array([-1.42108547e-14])
1
func(-6)
1
6
1
func(-7)
1
-70
1
func(-8)
1
-184
1
2
plt.plot(np.linspace(-9,5,20),func(np.linspace(-9,5,20)))
plt.axhline(0)
1
<matplotlib.lines.Line2D at 0x1fad4cb7700>

png

1
from scipy.optimize import least_squares
1

1

发表于 更新于
本文字数： 121 阅读时长 ≈ 1 分钟

markdown 测试文档

一级标题

二级标题

三级标题

四级标题

  • 无序列表
  • 无序列表
  • 无序列表
  1. 有序列表
  2. 有序列表

-[ ] 待办
-[x] 已完成

  • python 代码
1
print('hello,world')
  • 引用块

引用块

  • 公式

$\alpha = \beta$

发表于 更新于
本文字数： 1.1k 阅读时长 ≈ 1 分钟

数学建模

  1. 误差
  2. 插值函数
  • 多项式插值：
  • lagrange插值Lagrange多项式
  • 两点插值（线性插值）
  • 三点插值（二次插值；抛物线插值）
    不具有计算继承性：如果增减点，旧的计算就得重算
  • Newton插值Newton多项式

将求解N+1元方程组的问题转化为求解n+1个一元一次方程的问题
具有继承性
以上三种插值即选取了三组不同的基函数。

  • 等距节点的Newton插值公式
    • 向前差分:差分表
    • $a_k = \frac{\Delta^k(y_0)}{k!h^k}$
    • Newton 向前插值公式（前插公式）
    • 另有，向后差分，中心差分，
    • 不等距节点Newton插值公式—差商：差商表
    • 差商的一般性质
  • 插值多项式的误差
    • 插值余项

平面方程

  • 点法式方程
  • 一般方程
  • 截距式方程
  • 两平面的夹角
  • 点到直线的距离

设平面上一点$M_0(x_0,y_0,z_0)$,平面法向量为$\vec{\bm{n}}=(A,B,C)$;

点法式方程：

$A(x-x_0)+B(y-y_0)+C(z-z_0)=0$

一般方程：

$Ax+By+Cz+D=0$

  • D=0,平面通过坐标原点；
  • A=0,平面平行于x轴；
  • A=B=0,平面平行于xoy面(垂直于z轴)；
  • A=D=0,平面通过x轴。

截距式方程：

设平面与x,y,z轴分别交于点$(a,0,0),(0,b,0),(0,0,c)$，则

${x\over a}+{y\over b}+{z\over c}=1$

两平面的夹角：

$\cos \varphi = {|A_1 A_2+B_1 B_2+C_1 C_2|\over \sqrt{A_2^2+B_2^2+C_2^2}\cdot \sqrt{A_1^2+B_1^2+C_1^2}}$

  • 两平面垂直

$A_1 A_2+B_1 B_2+C_1 C_2= 0$

  • 两平面平行

${A_1\over A_2}= {B_1\over B_2}={C_1\over C_2}$

点到平面的距离：

$d = {|Ax_0+By_0+Cz_0+D|\over\sqrt{A^2+B^2+C^2}}$

发表于 更新于
本文字数： 367 阅读时长 ≈ 1 分钟

Welcome to Hexo! This is your very first post. Check documentation for more info. If you get any problems when using Hexo, you can find the answer in troubleshooting or you can ask me on GitHub.

Quick Start

Create a new post

1
$ hexo new "My New Post"

More info: Writing

Run server

1
$ hexo server

More info: Server

Generate static files

1
$ hexo generate

More info: Generating

Deploy to remote sites

1
$ hexo deploy

More info: Deployment

发表于 更新于
本文字数： 48k 阅读时长 ≈ 43 分钟

《科学计算与matlab语言》，中国大学MOOC， 中南大学

https://www.icourse163.org/learn/CSU-1002475002#/learn/content

使用中的技巧发现

Pwd 输出当前工作目录，支持Linux命令。

diary - 将命令行窗口文本记录到日志文件中

save - 将工作区变量保存到文件中

load - 将文件变量加载到工作区中

0 课程导入

求x^2-3x+1=0的根。

方法一：利用MATLAB多项式求根函数roots来求根。

1
2
p=[1,-3,1];
x=roots(p)

绘图：

1
2
3
4
x=-5:0.1:5;
y1=x.*x-3*x+1;
y2=zeros(size(x));
plot(x, y1, x, y2);

方法二 ： 利用求单变量非线性方程根的函数fzero，求方程在某个初始点附近的实根。

1
2
3
4
 f=@(x) x*x-3*x+1;
 x1=fzero(f, 0.5)

x2=fzero(f, 2.5)

方法三：利用最优化工具箱中的方程求根函数fsolve。

1
2
3
4
 f=@(x) x*x-3*x+1;
 x1=fsolve(f, 0.5, optimset('Display', 'off'))

x2=fsolve(f, 2.5, optimset('Display', 'off'))

方法四：利用solve函数求方程的符号解，即求得的解是一个表达式。

1
2
3
4
 syms x
 x=solve(x^2-3*x+1)

x=eval(x)

1 matlab基础知识

1.1 matlab系统环境

1.2 matlab数值数据

例1 分别求一个三位正整数的个位数字、十位数字和百位数字。

1
2
3
4
m=345;
m1=rem(m,10)
m2=rem(fix(m/10),10)
m3=fix(m/100)

例2 求[1，100]区间的所有素数。

1
2
3
4
x=1:100;
k=isprime(x);
k1=find(k);
p=x(k1)

1.3 变量及其操作

例1 计算表达式的值，并将结果赋给变量z，然后显示计算结果。

1
2
3
x=sqrt(7)-2i;
y=exp(pi/2);
z=(5+cosd(47))/(1+abs(x-y))

1.4 MATLAB矩阵的表示方法

1.5 矩阵元素的引用

1.6 MATLAB基本运算

例1 当x=0.1、0.4、0.7、1时，分别求y=sin x cos x的值。

1
2
x=0.1:0.3:1;
y=sin(x).*cos(x)

例2 建立3阶方阵A，判断A的元素是否为偶数。

1
2
A =[24,35,13;22,63,23;39,47,80]
P=rem(A,2)==0

例3 水仙花数是指各位数字的立方之和等于该数本身的三位正整数。求全部水仙花数。

1
2
3
4
5
6
m=100:999;
m1=rem(m,10);
m2=rem(fix(m/10),10);
m3=fix(m/100);
k=find(m==m1.*m1.*m1+m2.*m2.*m2+m3.*m3.*m3)
s=m(k)

1.7 字符串处理

例1 建立一个字符串向量，然后对该向量做如下处理：
① 取第1～5个字符组成的子字符串。
② 将字符串倒过来重新排列。
③ 将字符串中的小写字母变成相应的大写字母，其余字符不变。
④ 统计字符串中小写字母的个数。

1
2
3
4
5
6
ch='ABc123d4e56Fg9';
subch=ch(1:5)  
revch=ch(end:-1:1) 
k=find(ch>='a'&ch<='z') 
ch(k)=ch(k)-('a'-'A')    
length(k)

专题一总结

2 matlab 矩阵处理

2.1 特殊矩阵

例1 产生5阶两位随机整数矩阵A，再产生均值为0.6、方差为0.1的5阶正态分布随机矩阵B，验证(A+B)I=IA+BI（I为单位矩阵）。

1
2
3
4
A=fix(10+(99-10+1)*rand(5));
B=0.6+sqrt(0.1)*randn(5);
C=eye(5);
(A+B)*C==C*A+B*C

例2 产生8阶魔方阵，求其每行每列元素的和。

1
2
3
M=magic(8);
sum(M(1,:))
sum(M(:,1))

例3 生成5阶帕斯卡矩阵，验证它的逆矩阵的所有元素也为整数。

1
2
3
format rat
P=pascal(5)
inv(P)

2.2 矩阵变换

例1 先建立5×5矩阵A，然后将A的第一行元素乘以1，第二行乘以2，…，第五行乘以5。

1
2
3
A=[7,0,1,0,5;3,5,7,4,1;4,0,3,0,2;1,1,9,2,3;1,8,5,2,9]
D=diag(1:5);
D*A

要将A的各列元素分别乘以对角阵的对角线元素。

1
2
3
A=[7,0,1,0,5;3,5,7,4,1;4,0,3,0,2;1,1,9,2,3;1,8,5,2,9]
D=diag(1:5);
A*D

例2 验证魔方阵的主对角线、副对角线元素之和相等。

1
2
3
4
5
6
7
A=magic(5);
D1=diag(A);
sum(D1)

B=flipud(A);
D2=diag(B);
sum(D2)

例3 用求逆矩阵的方法解线性方程组Ax=b。

1
2
3
4
A=[1,2,3;1,4,9;1,8,27];
b=[5;-2;6];
x=inv(A)*b
x=A\b

2.3 矩阵求值

例1 验证det(A-1)=1/det(A)。

1
2
3
4
format rat
A=[1,3,2;-3,2,1;4,1,2]
det(inv(A))
1/det(A)

例2 求3~20阶魔方阵的秩。

1
2
3
4
5
6
7
for n=3:20
  r(n)=rank(magic(n));
end
bar(r)
grid on
axis([2,21,0,20])
[3:20;r(3:20)]

例3 求2~10阶希尔伯特矩阵的条件数。

1
2
3
4
5
for n=2:10
  c(n)=cond(hilb(n));
end
format long
c'

2.4 矩阵的特征值与特征向量

例1 

1
2
3
4
5
6
R=[-1,2,0;2,-4,1;1,1,-6];
S=[1,2;2,3];
A=[R,zeros(3,2);zeros(2,3),S];
[X1,d1]=eig(R)
[X2,d2]=eig(S)
[X3,d3]=eig(A)

例2

1
2
3
4
5
6
7
x=[0,0.5,0.5,3,5.5,5.5,6,6,3,0;0,0,6,0,6,0,0,8,1,8];
A=[1,0.5;0,1];
y=A*x;
subplot(2,2,1);
fill(x(1,:),x(2,:),'r');
subplot(2,2,2);
fill(y(1,:),y(2,:),'r');

2.5 稀疏矩阵

求下列三对角线性方程组的解。

1
2
3
4
5
6
7
8
kf1=[1;1;2;1;0];
k0=[2;4;6;6;1];
k1=[0;3;1;4;2];
B=[kf1,k0,k1];
d=[-1;0;1];
A=spdiags(B,d,5,5);
f=[0;3;2;1;5];
x=A\f

专题二总结

3 matlab 程序流程控制

3.1 顺序结构程序

例1 有一线段AB，A的坐标为(1，1)，B的坐标为(4.5，4.5)，求AB的长度，以及黄金分割点C的坐标。

1
2
3
4
5
6
a=input('a=');
b=input('b=');
c=a+0.618*(b-a);
s=abs(a-b);
disp(s)
disp(c)

3.2 用if语句实现选择结构

例1 输入一个整数，若为奇数则输出其平方根，否则输出其立方根。 

1
2
3
4
5
6
7
x=input('请输入x的值:');
if rem(x,2)==1
  y=sqrt(x);
else
  y=x^(1/3);
end
y

例2 输入一个字符，若为大写字母，则输出其对应的小写字母；若为小写字母，则输出其对应的大写字母；若为数字字符则输出其对应数的平方，若为其他字符则原样输出。

1
2
3
4
5
6
7
8
9
10
c=input('请输入一个字符：','s');
if c>='A' && c<='Z'
  disp(lower(c))
elseif c>='a' && c<='z'
  disp(upper(c))
elseif c>='0' && c<='9'
  disp(str2double(c)^2)
else
  disp(c)
end

3.3 用switch语句实现选择结构

例1 输入一个英文单词，判断它是否以元音字母开头。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
c=input('请输入一个单词：','s');
switch c(1) 
  case {'A','E','I','O','U','a','e','i','o','u'} 
    disp([c,'以元音字母开头']);
  otherwise 
    disp([c,'以辅音字母开头']);
end

c=input('请输入一个单词：','s');
if findstr(c(1),'AEIOUaeiou')>0
  disp([c,'以元音字母开头']);
else
  disp([c,'以辅音字母开头']);
end

例2 PM2.5是指大气中直径小于或等于2.5微米的可入肺颗粒物，是衡量空气质量的重要指标。假定空气质量等级以PM2.5数值划分为6级。 PM2.5数值[0，35）空气质量为优，[35，75）为良，[75，115）为轻度污染，[115，150）为中度污染，[150，250）为重度污染，大于等于250为严重污染。编写程序，输入PM2.5数值，输出空气质量等级。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
g=input('请输入PM2.5值：');
switch fix(g) 
  case num2cell(0:34) 
    disp('空气质量优');
  case num2cell(35:74) 
    disp('空气质量良好');
  case num2cell(75:114) 
    disp('空气质量轻度污染');
  case num2cell(115:149) 
    disp('空气质量中度污染'); 
  case num2cell(150:249) 
    disp('空气质量重度污染'); 
  otherwise 
    disp('空气质量严重污染');
end

3.4 用for语句实现循环

计算圆周率π。

（1）利用无穷级数展开式求π的近似值。

1
2
3
4
5
6
7
8
y=0;
g=-1;
n=input('n=?'); 
for i=1:n
  g=-g;
  y=y+g/(2*i-1); 
end
pai=4*y

用向量求和的方法实现程序：

1
2
3
4
n=input('n=?');
x=1:2:(2*n-1);
y=(-1).^(2:n+1)./x;
pai=sum(y)*4

（2）利用定积分的近似值求π的近似值。

1
2
3
4
5
6
7
8
9
10
11
12
a=0;
b=1;
n=input('n=?');
h=(b-a)/n;
x=a:h:b;
f=sqrt(1-x.*x);
s=[];
for k=1:n
  s1=(f(k)+f(k+1))*h/2;
  s=[s,s1];
end
pai=4*sum(s)

（3）利用蒙特卡洛法求π的近似值。

1
2
3
4
5
6
7
8
9
10
11
12
>  s=0;
>  n=input('n=?');
>  for i=1:n
>    x=rand(1);
>    y=rand(1);
>    if x*x+y*y<=1
>      s=s+1;
>    end
>  end
>  pai=s/n*4

3.5 用while语句实现循环

例1 从键盘输入若干个数，当输入0时结束输入，求这些数的平均值和它们之和。

1
2
3
4
5
6
7
8
9
10
11
12
msum=0;
n=0;
x=input('Enter a number (end in 0):');
while x~=0
  msum=msum+x;
  n=n+1;
  x=input('Enter a number (end in 0):');
end
if n>0
  msum
  mean=msum/n
end

例2 求[100，200]之间第一个能被21整除的整数。

1
2
3
4
5
6
7
for n=100:200
  if rem(n,21)~=0
    continue
  end
  n
  break
end

例3 用筛选法求某自然数范围内的全部素数。

1
2
3
4
5
6
7
8
9
10
m=input('m='); 
p=1:m;
p(1)=0;
for i=2:sqrt(m)
  for j=2*i:i:m
    p(j)=0;
  end
end
n=find(p~=0);
p(n)

3.6 函数文件的定义与调用

例1 编写函数文件，求半径为r的圆的面积和周长。

1
2
3
function [s,p]=fcircle(r)
s=pi*r*r;
p=2*pi*r;

例2 已知=(40)/((30)+(20))
①当()=+10 ln⁡(^2+5)时，的值是多少。
②当()=1×2+2×3+3×4+⋯+×(+1)时，的值是多少。
分别用匿名函数和函数文件定义函数()。
第②问的函数文件f2.m。

1
2
3
4
5
function f=f2(n)
f=0;
for k=1:n
  f=f+k*(k+1);
end

脚本文件mf.m。

1
2
3
f1=@(n) n+10*log(n*n+5);
y1=f1(40)/(f1(30)+f1(20))
y2=f2(40)/(f2(30)+f2(20))

3.7 函数的递归调用

例1 利用函数的递归调用，求n！。

函数文件fact.m如下：

1
2
3
4
5
6
function f=fact (n)
if n<=1
  f=1;
else
  f=fact (n-1)*n;  %递归调用求(n-1)!
end

在脚本文件a.m中调用函数文件fact.m，求n!。

1
2
3
n=input('Please input n=');
s=fact (n);
disp(s)

例2 Fibonacci数列定义如下：
f1=1
f2=1
fn=fn-1+fn-2 （n>2）
编写递归调用函数求Fibonacci数列的第n项，然后调用该函数验证Fibonacci数列的如下性质：

首先建立函数文件ffib.m。

1
2
3
4
5
6
function f=ffib(n)
if n>2
  f=ffib(n-1)+ffib(n-2);
else
  f=1;
end

建立程序文件test.m。

1
2
3
4
5
6
F=[];
for k=1:20
  F=[F,ffib(k)*ffib(k)];
end
sum(F)
ffib(20)*ffib(21)

3.8 函数参数与变量作用范围

例1 利用nargin和nargout建立函数文件test.m。

1
2
3
4
5
6
7
8
function fout=test(a,b,c)
if nargin==1
  fout=a;
elseif nargin==2
  fout=a+b;
elseif nargin==3
  fout=(a*b*c)/2;
end

例2 利用全局变量建立函数文件wad.m。

1
2
3
function f=wad(x,y)
global ALPHA BETA
f=ALPHA*x+BETA*y;

专题三总结

4 matlab 绘图

4.1 二维曲线

例1 绘制一条折线。

1
2
3
x=[2.5, 3.5, 4, 5];
y=[1.5, 2.0, 1, 1.5];
plot(x, y)

例2 绘制sin(x)、sin(2x)、sin(x/2)的函数曲线。

1
2
3
4
5
6
7
8
9
 x=linspace(0,2*pi,100);
 y=[sin(x); sin(2*x); sin(0.5*x)];
 plot(x,y)

t=0:0.01:2*pi;
 t1=t';
 x=[t1, t1, t1];
 y=[sin(t1), sin(2*t1), sin(0.5*t1)];
 plot(x,y)

例3 采用不同个数的数据点绘制正弦函数曲线，观察曲线形态。

1
2
3
4
t1=linspace(0, 2*pi, 10);
t2=linspace(0, 2*pi, 20);
t3=linspace(0, 2*pi, 100);
plot(t1, sin(t1), t2,sin(t2)+1, t3, sin(t3)+2)

例4 用不同线型和颜色在同一坐标内绘制曲线及其包络线。

1
2
3
4
5
6
x=(0:pi/50:2*pi)';
y1=2*exp(-0.5*x)*[1,-1];
y2=2*exp(-0.5*x).*sin(2*pi*x);
x1=0:0.5:6;
y3=2*exp(-0.5*x1).*sin(2*pi*x1);
plot(x, y1, 'k:', x, y2, 'b--', x1, y3, 'rp')

例5 绘制函数sin(1/)的图形。

1
2
3
x=0:0.005:0.2;
y=sin(1./x);
plot(x,y)

例6 采用fplot函数绘制函数sin⁡(1/x)。

1
fplot(@(x) sin(1./x),[0,0.2], 'b')

例7 已知螺旋线的参数方程，绘制曲线。

1
fplot(@(t)t.*sin(t), @(t)t.*cos(t), [0,10*pi], 'r')

4.2 绘制图形的辅助标注

例1 绘制[-2π，2π]区间的正弦曲线并给图形添加标题。

1
2
3
4
5
x=-2*pi:0.05:2*pi;
y=sin(x);
plot(x,y)
title('y=sin(x)')
title({'MATLAB', 'y=sin(x)'})

\bf 加粗

\it 斜体

\rm 正体

Color（’k’）,fontsize(11)

绘制[-2π，2π]区间的正弦曲线并给x轴添加标签。

1
2
3
4
5
x=-2*pi:0.05:2*pi;
y=sin(x);
plot(x,y)
title('y=sin(x)')
xlabel('-2\pi \leq x \leq 2\pi')

在前面的图形中添加文字说明。Text(x,y,txt);gtext(txt)

1
2
text(-2*pi, 0, '-2{\pi}')
 text(3, 0.28, '\leftarrow sin(x)')

例2 绘制不同频率的正弦曲线并用图例标注曲线。Legend()

1
2
3
x = linspace(0, 2*pi, 100);
plot(x, [sin(x); sin(2*x); sin(3*x)])
legend('sin(x)', 'sin(2x)', 'sin(3x)')

绘制边长为1的正方形，并设置坐标轴。Axis()

1
2
3
4
5
x = [0, 1, 1, 0, 0];
 y = [0, 0, 1, 1, 0];
 plot(x,y)
 axis([-0.1, 1.1, -0.1, 1.1])
 axis equal;

axis equal

axis square

axis auto

axis off

axis on

grid on

grid off

grid

box on

box off

box

例3 绘制sin(x)、sin(2x)、sin(x/2)的函数曲线并添加图形标注。

1
2
3
4
5
6
7
8
9
10
11
x=linspace(0,2*pi,100);
y=[sin(x); sin(2*x); sin(0.5*x)];
plot(x,y)
axis([0 7 -1.2, 1.2]) 
title('不同频率正弦函数曲线'); 
xlabel('Variable X'); ylabel('Variable Y'); 
text(2.5, sin(2.5), 'sin(x)'); 
text(1.5, sin(2*1.5), 'sin(2x)');
text(5.5, sin(0.5*5.5), 'sin(0.5x)');
legend('sin(x)','sin(2x)','sin(0.5x)')
grid on

例4 用图形保持功能绘制两个同心圆。Hold on ;hold off;hold

1
2
3
4
5
6
7
8
t = linspace(0,2*pi,100);
x = sin(t); y = cos(t);
plot(x, y, 'b')
hold on;    
plot(2*x, 2*y, 'r--')
grid on
axis([-2.2 2.2 -2.2 2.2])
axis equal

划分2×2子图 subplot(m,n,p)

1
2
3
4
5
6
subplot(2,2,1); 
x=linspace(0,2*pi,60);
y=sin(x);
plot(x,y);
title('sin(x)');
axis ([0,2*pi,-1,1]);

划分多子图

1
2
3
4
5
6
7
8
9
10
11
12
13
x=linspace(0,2*pi,60);
subplot(2,2,1)   
plot(x,sin(x)-1);
title('sin(x)-1');axis ([0,2*pi,-2,0])
subplot(2,1,2)   
plot(x,cos(x)+1);
title('cos(x)+1');axis ([0,2*pi,0,2])
subplot(4,4,3)   
plot(x,tan(x));
title('tan(x)');axis ([0,2*pi,-40,40])
subplot(4,4,8)    
plot(x,cot(x));
title('cot(x)');axis ([0,2*pi,-35,35])

4.3 其他形式的二维图形

例1 绘制1/的直角线性坐标图和三种对数坐标图。

对数坐标图semilogx();seilogy();loglog()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
x=0:0.1:10;
y=1./x;
subplot(2,2,1)
plot(x,y) 
title('plot(x,y)');grid on
subplot(2,2,2)
semilogx(x,y)  
title('semilogx(x,y)');grid on
subplot(2,2,3)
semilogy(x,y) 
title('semilogy(x,y)');grid on
subplot(2,2,4)
loglog(x,y) 
title('loglog(x,y)');grid on

例2 按极坐标方程ρ=1-sin t绘制心形曲线。Polar(theta,rho,options)

1
2
3
4
5
6
7
8
t = 0:pi/100:2*pi;
r = 1-sin(t);
subplot(1,2,1)
polar(t,r)
subplot(1,2,2)
t1 = t-pi/2;
r1 = 1-sin(t1);
polar(t,r1)

例3 绘制分组条形图。bar(y,style);barh();stacked:堆积分组；groud:簇状分组（默认）

1
2
3
4
5
6
7
y=[1,2,3,4,5; 1,2,1,2,1; 5,4,3,2,1];
subplot(1,2,1)
bar(y)
title('Group')
subplot(1,2,2)
bar(y, 'stacked')
title('Stack')

例4 下表是某公司2015~2017年家电类商品1月份的销售数据，绘制条形图对比数据。

1
2
3
4
5
6
x=[2015,2016,2017];
y=[68,80,115,98,102;
75,88,102,99,110;
81,86,125,105,115];
bar(x, y)
title('Group');

例5 绘制服从高斯分布的直方图。hist();直角坐标系；rose():极坐标系。

1
2
3
4
5
6
7
8
y=randn(500,1);
subplot(2,1,1); 
hist(y); 
title('高斯分布直方图');
subplot(2,1,2); 
x=-3:0.2:3;
hist(y,x); 
title('指定区间中心点的直方图')')

例6 绘制高斯分布数据在极坐标下的直方图。

1
2
3
4
y=randn(500,1);
theta=y*pi;
rose(theta)
title('在极坐标下的直方图')

例7 某次考试优秀、良好、中等、及格、不及格的人数分别为：7、17、23、9、4，试用扇形统计图作成绩统计分析。Pie函数：扇形图；area函数：面积图

1
2
3
4
5
score = [5, 17, 23, 9, 4];
ex = [0,0,0,0,1];
pie(score, ex)  
legend('优秀', '良好', '中等’, '及格', '不及格', …
'location', 'eastoutside')

例8 以散点图形式绘制桃心曲线。Scatter散点图；stairs阶梯图;stem杆图

1
2
3
4
t = 0:pi/50:2*pi;
x = 16*sin(t).^3;
y = 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t);
scatter(x,y,'rd','filled')

例9 已知向量A、B，求A+B，并用矢量图表示。Compass 罗盘图；feather羽毛图；quiver箭头图

1
2
3
4
5
6
7
8
A=[4,5]; B=[-10,0]; C=A+B;
hold on;
quiver(0, 0, A(1), A(2));
quiver(0, 0, B(1), B(2));
quiver(0, 0, C(1), C(2));
text(A(1),A(2),'A');text(B(1),B(2),'B'); text(C(1),C(2),'C'); 
axis ([-12, 6, -1, 6])
grid on

4.4 三维曲线

例1 绘制一条空间折线。Plot3(x,y,z);fplot3

1
2
3
4
5
6
x=[0.2, 1.8, 2.5];
y=[1.3, 2.8, 1.1];
z=[0.4, 1.2, 1.6];
plot3(x, y, z)
grid on
axis([0, 3, 1, 3, 0, 2]);

例2 绘制螺旋线 

1
2
3
4
5
6
7
8
9
10
t=linspace(0, 10*pi, 200);
x=sin(t)+t.*cos(t);
y=cos(t)-t.*sin(t);
z=t;
subplot(1, 2, 1)
plot3(x, y, z)
grid on
subplot(1, 2, 2)
plot3(x(1:4:200), y(1:4:200), z(1:4:200))
grid on

例3 在空间不同位置绘制5条正弦曲线。

1
2
3
4
5
6
7
8
9
10
11
12
t=0:0.01:2*pi;
t=t';
x=[t, t, t, t, t];
y=[sin(t), sin(t)+1, sin(t)+2, sin(t)+3, sin(t)+4];
z=x;
plot3(x,y,z)
这个例子也可以采用以下代码实现。
t=0:0.01:2*pi;
x=t;
y=[sin(t); sin(t)+1; sin(t)+2; sin(t)+3; sin(t)+4];
z=x;
plot3(x,y,z)

例4 绘制三条不同长度的正弦曲线。 

1
2
3
4
5
t1=0:0.01:1.5*pi;
t2=0:0.01:2*pi;
t3=0:0.01:3*pi;
plot3(t1,sin(t1),t1, t2,sin(t2)+1,t2, …
t3,sin(t3)+2,t3)

例5 绘制空间曲线

1
2
3
4
5
6
7
t=0:pi/50:6*pi;
x=cos(t);
y=sin(t);
z=2*t;
plot3(x,y,z,'p')
xlabel('X'),ylabel('Y'),zlabel('Z');
grid on

例6 绘制墨西哥帽顶曲线。参数方程fplot3(funx,funy,funz,tlims[-5,5])

1
2
3
4
xt = @(t) exp(-t/10).*sin(5*t);
yt = @(t) exp(-t/10).*cos(5*t);
zt = @(t) t;
fplot3(xt, yt, zt, [-12, 12])

用红色点划线绘制墨西哥帽顶曲线。

1
2
3
4
xt = @(t) exp(-t/10).*sin(5*t);
yt = @(t) exp(-t/10).*cos(5*t);
zt = @(t) t;
fplot3(xt, yt, zt, [-12, 12], 'r-.')

4.5 三维曲面：

mesh三维网格图;surf三维曲面图;fmesh;fsurf;;meshgrid；meshc;meshz;surfc;surfl;

例1 绘制空间曲线。

1
2
3
4
5
6
x = 2:6; 
y = (3:8)';
[X, Y] = meshgrid(x, y);
Z = randn(size(X));
plot3(X,Y,Z)
grid on;

例2 绘制三维曲面图。 

1
2
3
4
5
6
7
8
9
10
t = -2:0.2:2; 
[X, Y] = meshgrid(t);
Z = X .* exp(-X.^2 - Y.^2);
subplot(1,3,1)
mesh(X,Y,Z); 
subplot(1,3,2)
surf(X,Y,Z); 
subplot(1,3,3)
plot3(X,Y,Z); 
grid on

例3 用4种方式绘制函数$z=(x-1)^2+(y-2)^2-1$曲面图。
其中，x∈[0,2]，y∈[1,3]。

1
2
3
4
5
6
7
8
9
10
[x,y]=meshgrid(0:0.1:2,1:0.1:3);
 z=(x-1).^2+(y-2).^2-1;
 subplot(2,2,1);
 meshc(x,y,z);title('meshc(x,y,z)')
 subplot(2,2,2);
 meshz(x,y,z);title('meshz(x,y,z)')
 subplot(2,2,3);
 surfc(x,y,z);title('surfc(x,y,z)')
 subplot(2,2,4);
 surfl(x,y,z); title('surfl(x,y,z)')

例4 用cylinder函数分别绘制柱面、花瓶和圆锥面。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Sphere(n=20);cylinder(R,n=20)
 subplot(1,3,1);
 [x,y,z]=cylinder;
 surf(x,y,z);
 subplot(1,3,2);
 t=linspace(0,2*pi,40);
 [x,y,z]= cylinder(2+cos(t),30);
 surf(x,y,z);
 subplot(1,3,3);
 [x,y,z]= cylinder(0:0.2:2,30);
 surf(x,y,z);

 

Peaks()

例5 用cylinder函数绘制两个相互垂直且直径相等的圆柱面的相交图形。

1
2
3
4
5
6
[x,y,z]= cylinder(1,60);
z=[-1*z(2,:);z(2,:)];
surf(x,y,z)
hold on
surf(y,z,x)
axis equal

例6 绘制螺旋曲面。

1
2
3
4
5
6
7
funx = @(u,v) u.*sin(v);
funy = @(u,v) -u.*cos(v);
funz = @(u,v) v;
fsurf(funx,funy,funz,[-5 5 -5 -2])
hold on
fmesh(funx,funy,funz,[-5 5 -2 2])
hold off

4.6 图形修饰处理

视点处理：方位角；仰角；view(az=-37.5,el=30)

View(x,y,z)

色彩处理[RGB]

colormap

剪裁处理

例1 绘制函数$z=(x-1)^2+(y-2)^2-1$曲面，并从不同视点展示曲面。

1
2
3
4
5
6
7
8
9
10
[x,y]=meshgrid(0:0.1:2, 1:0.1:3);
z=(x-1).^2+(y-2).^2-1;
subplot(2,2,1); mesh(x,y,z)
title('方位角=-37.5{\circ},仰角=30{\circ}')
subplot(2,2,2); mesh(x,y,z)
view(0,90);title('方位角=0{\circ},仰角=90{\circ}')
subplot(2,2,3); mesh(x,y,z)
view(90,0); title('方位角=90{\circ},仰角=0{\circ}')
subplot(2,2,4); mesh(x,y,z)
view(-45,-60); title('方位角=-45{\circ},仰角=-60{\circ}')

例2 创建一个灰色系列色图矩阵。

1
2
3
4
5
6
7
8
c = [0, 0.2, 0.4, 0.6, 0.8, 1.0]';
cmap = [c, c, c];
surf(peaks)
colormap(cmap)
或
cmap=gray(6);
surf(peaks)
colormap(cmap)

例3 使用同一色图，以不同着色方式绘制圆锥体。Shading faceted;shading flat;shading interp

1
2
3
4
5
6
7
8
[x,y,z]= cylinder(pi:-pi/5:0,10);
colormap(lines);
subplot(1,3,1);
surf(x,y,z); shading flat
subplot(1,3,2);
surf(x,y,z); shading interp
subplot(1,3,3);
surf(x,y,z);

例4 绘制3/4圆。将图形中需要裁减部分对应的函数值设置成NaN。

1
2
3
4
5
6
7
8
9
t = linspace(0,2*pi,100);
x = sin(t);
y = cos(t);
p = y > 0.5;
y(p)= NaN;
plot(x,y)
axis([-1.1,1.1,-1.1,1.1])
axis square
grid on

例5 绘制3/4球面。

1
2
3
4
5
6
7
[X, Y, Z] = sphere(60);
p = Z>0.5;
Z(p) = NaN;
surf(X, Y, Z)
axis([-1, 1, -1, 1, -1, 1])
axis equal
view(-45, 20)

4.7 交互式绘图工具

“绘图”选项卡

图形窗口绘图工具：plottools;按钮；；图形选项板（新子图面板、变量面板、注释面板）、绘图浏览器、属性编辑器

图形窗口菜单和工具栏（工具栏、图形窗口菜单）

生成例1所需变量。

1
2
3
4
5
6
7
8
x=linspace(0,2*pi,100);
y=sin(x);
y1=sin(x);
y2=sin(0.5*x);
y3=sin(2*x);
[u,v]=meshgrid(0:0.1:2,1:0.1:3);
h=(u-1).^2+(v-2).^2-1;
[X,Y,Z]= cylinder(0:0.2:2,30);

专题四总结

5、数据分析与多项式计算

5.1 数据统计分析(matlab以列优先)

Max/min/ 最大值、最小值

Mean/median 均值、中值

Sum/prod 求和、求积

Cumsum/cumprod 累加和、累乘积

Std/corrcoef 标准差、相关系数

样本标准差、总体标准差

Sort 排序

例1 求向量x的最大元素，其中x=[-43,72,9,16,23,47]。

1
2
3
4
\>> x=[-43,72,9,16,23,47];
\>> y=max(x) 

\>> [y,k]=max(x)

例2 求矩阵A（见课件）的每行及每列的最大元素，并求整个矩阵的最大元素。

1
2
3
4
5
6
\>> A=[13,-56,78;25,63,-235;78,25,563;1,0,-1];
\>> max(A)

\>> max(A,[],2)

\>> max(max(A))

例3 某学生宿舍的5位同学月生活费如向量x所示，其中，小明同学家境一般，请问他应该按什么标准向父母主张生活费额度才较为合理。x=[1200,800,1500,1000,5000]

1
2
3
4
\>> x=[1200,800,1500,1000,5000];
\>> mean(x)

\>> median(x)

例4 求向量X=[1,2,3,4,5,6,7,8,9,10]的积与累乘积。

1
2
3
4
\>> X=[1,2,3,4,5,6,7,8,9,10];
\>> y1=prod(X)

\>> y2=cumprod(X)

例5 生成满足正态分布的50000*4随机矩阵，用不同的形式求其各列之间的标准差。

1
2
3
4
5
6
7
8
9
10
11
\>> x=randn(50000,4);
\>> y1=std(x,0,1)

\>> y2=std(x,1,1)

\>> x1=x';
\>> y3=std(x1,0,2);
\>> y3'

\>> y4=std(x1,1,2);
\>> y4'

例6 某新产品上市，在上市之前，公司物流部门把新产品分配到不同地区的10个仓库进行销售。产品上市一个月后，公司要对各种不同的分配方案进行评估，以便在下一次新产品上市时进行更准确的分配，避免由于分配不当而产生的积压和断货。下表（见课件）是相关数据，请判断那种分配方案最为合理。

1
2
3
4
5
6
\>> A=[5032,6000,5100,5200;6532,6500,6600,5800;
   5500,7000,5400,4800;4530,4000,4300,4200;
   2300,2000,2200,2500;3254,3000,3500,3000;
   8095,9000,7800,8500;7530,8000,7000,7500;
   3841,3200,3500,3200;4500,5200,4800,4000];
\>> corrcoef(A)

例7 对下列矩阵（见课件）做各种排序。

1
2
3
4
5
6
\>> A=[1,-8,5;4,12,6;13,7,-13];
\>> sort(A)

\>> sort(A,2,'descend')

\>> [X,I]=sort(A)

5.2 多项式计算

多项式的表示：n+1长度的行向量

多项式的四则预算：加减（相应向量加减）

P=Conv(p1,p2);相乘

[Q,r]=Deconv(p1,p2);除法;p1 = conv(Q,p2)+r

多项式的求导：polyder;

多项式的求值：polyval(p,x);ployvalm(p,x)

矩阵运算、点运算

多项式的求根：roots(p);给出全部根；poly函数由全部根建立多项式

例1 设f和g为两个多项式（见课件），求f(x)+g(x),f(x)-g(x),f(x)×g(x),f(x)/g(x)。

1
2
3
4
5
6
7
8
9
10
\>> f=[3,-5,0,-7,5,6]; g=[3,5,-3]; g1=[0,0,0,g];
\>> f+g1

\>> f-g1

\>> conv(f,g)

\>> [Q,r]=deconv(f,g)

\>> conv(g,Q)+r

例2 已知两个多项式a和b（见课件），计算两个多项式的乘积的导函数、商的导函数。

1
2
3
4
5
6
7
 \>> a=[3 1 0 -6];
 \>> b=[1 2];
 \>> polyder(a)

\>> c=polyder(a,b)
 
 \>> [p,q]=polyder(a,b)

>> [p,q]=polyder(a,b)

例3 以一个多项式为例，取一个2×2矩阵为自变量，分别用polyval和polyvalm计算该多项式的值。

1
2
3
4
5
\>> a=[1,8,0,0,-10];
\>> x=[-1,1.2;2,-1.8];
\>> y1=polyval(a,x)

\>> y2=polyvalm(a,x)

例4 请推算空气流量在[0, 2]范围内什么水平时，加热效率为最高。

1
2
3
4
5
6
7
8
9
\>> p=[-38.89,126.11,-3.42];
\>> q=polyder(p)

\>> roots(q)

\>> polyval(p,1.6214)

\>> x=0:0.1:2;
\>> plot(x,polyval(p,x),1.6214,98.8154,'rp');

5.3 数据插值(函数逼近的方法)

数值计算方法

Interp1(x,y,x1,method)

Method:linear（线性插值）;nearest（最近点插值）;pchip（分段3次埃尔米特插值）;spline（3次样条插值）;

Interp2(x,y,z,x1,y1,method)

引例 零件加工问题 

1
2
3
4
5
\>> x= [0,3,5,7,9,11,12,13,14,15];
\>> y=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6];
\>> x1=0:0.1:15;
\>> y1=interp1(x,y,x1,'spline');
\>> plot(x1,y1)

粮储仓的通风控制问题

1
2
3
4
5
6
7
8
>> x=20:10:90;

\>> y=(0:5:20)';
 \>> z=[8.9,10.32,11.3,12.5,13.9,15.3,17.8,21.3;8.7,10.8,11,12.1,13.2,14.8,16.55,20.8;8.3,9.65,10.88,12,13.2,14.6,16.4,20.5;8.1,9.4,10.7,11.9,13.1,14.5,16.2,20.3;8.1,9.2,10.8,12,13.2,14.8,16.9,20.9];
 \>> xi=20:90;
 \>> yi=(0:20)';
 \>> zi=interp2(x,y,z,xi,yi,'spline');
 \>> surf(xi,yi,zi)

5.4 数据插值应用举例

机动车刹车距离问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 v=20:10:150;
 vs=v.*(1000/3600);
 d1=10.*vs;
 d2=[3.15,7.08,12.59,19.68,28.34,38.57,50.4,63.75,78.71,95.22,113.29,132.93,154.12,176.87];
 d3=10;
 d=d1+d2+d3;
 vi=20:1:150;
 di=interp1(v,d,vi,'spline');
 x=abs(di-120);
 [y,i]=sort(x);
 vi(i(1))
 plot(vi,di,vi(i(1)),di(i(1)),'rp')

\>> j=find(vi==125);
 \>> di(j)
 \>> plot(vi,di,125,480.14,'rp')

沙盘制作问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
x=0:200:1800;
y=x';
z=[2000,2000,2001,1992,1954,1938,1972,1995,1999,1999;2000,2002,2006,1908,1533,1381,1728,1959,1998,2000;
 2000,2005,2043,1921,977,897,1310,1930,2003,2000;1997,1978,2009,2463,2374,1445,1931,2209,2050,2003;
 1992,1892,1566,1971,2768,2111,2653,2610,2121,2007;1991,1875,1511,1556,2221,1986,2660,2601,2119,2007;
 1996,1950,1797,2057,2849,2798,2608,2303,2052,2003;1999,1999,2079,2685,3390,3384,2781,2165,2016,2000;
 2000,2002,2043,2271,2668,2668,2277,2049,2003,2000;2000,2000,2004,2027,2067,2067,2027,2004,2000,2000];
surf(x,y,z);
x1=0:100:1800;
y1=x1';
z1=interp2(x,y,z,x1,y1,'spline');
surf(x1,y1,z1);
x2=0:50:1800;
y2=x2';
z2=interp2(x1,y1,z1,x2,y2,'spline');
surf(x2,y2,z2);
contour(x2,y2,z2,12)

5.5 曲线拟合（一种函数逼近的方法）

Polyfit：求得最小二乘拟合多项式系数

原理：函数逼近；多项式函数；误差最小（最小二乘法/最小平方法）；

曲线拟合的实现方法：问题分析；

（若前面系数为0，则表明拟合次数太高，需要纠正）

引例—人口增长预测问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
 x=1790:10:2010;
 y=[3.9,5.3,7.2,9.6,12.9,17.1,23.2,31.4,38.6, 50.2,63.0,76.0,92.0,105.7,122.8,131.7,150.7,179.3,203.2,226.5,248.7,281.4,308.7];
 plot(x,y,'*');
 p=polyfit(x,y,3);
 polyval(p,2020)
 plot(x,y,'*',x,polyval(p,x));

\>> polyval(p,2016)
 
 \>> x=1950:10:2010;
 \>> y=[150.7,179.3,203.2,226.5,248.7,281.4,308.7];
 \>> p=polyfit(x,y,3)
 
 \>> p=polyfit(x,y,2);
 \>> plot(x,y,'*',x,polyval(p,x))
 \>> polyval(p,2016)
 
 \>> polyval(p,2020)

>> x=1950:10:2010;
>> y=[150.7,179.3,203.2,226.5,248.7,281.4,308.7];
>> p=polyfit(x,y,3)

>> p=polyfit(x,y,2);
>> plot(x,y,’*’,x,polyval(p,x))
>> polyval(p,2016)

>> polyval(p,2020)

家庭储蓄规律问题

1
2
3
4
5
\>> x=[0.6,1.0,1.4,1.8,2.2,2.6,3.0,3.4,3.8,4 ];
\>> y=[0.08,0.22,0.31,0.4,0.48,0.56,0.67,0.75,0.8,1.0 ];
\>> p=polyfit(x,y,1)

\>> plot(x,y,'*',x,polyval(p,x))

5.6 曲线拟合应用举例

股票预测问题

1
2
3
4
5
6
7
8
9
10
\>> x=[2,3,4,5,8,9,10,11,12,15,16,17,18,19,22,23,24,25,26,29,30];
\>> y=[7.74,7.84,7.82,7.78,7.91,7.97,7.9,7.76,7.9,8.04,8.06,8.11,8.08,8.13,8.03,8.01,8.06,8.0,8.3,8.41,8.28];
\>> plot(x,y, '*');
\>> p=polyfit(x,y,3);
\>> plot(x,y,'*',x,polyval(p,x));
\>> x1=[31,32,33];
\>> xi=[x,x1];
\>> plot(x,y,'*',xi,polyval(p,xi));
\>> y1=[8.27,8.17,9.54];
\>> plot(x,y,'*',xi,polyval(p,xi),x1,y1,'rp');

算法的参数优化问题

1
2
3
4
5
6
7
8
9
10
11
12
13
x=0.03:0.03:0.3;
y1=[0.01,0.01,0.02,0.03,0.06,0.07,0.13,0.17,0.25,0.37];
y2=[0.85,0.76,0.68,0.62,0.56,0.52,0.49,0.46,0.43,0.39];
plot(x,y1,'*',x,y2,'o');
p1=polyfit(x,y1,2);
p2=polyfit(x,y2,2);
p=p1-p2;
xi=roots(p);
xj=0:0.03:0.36;
yj1=polyval(p1,xj);
yj2=polyval(p2,xj);
yi=polyval(p1,xi(2))
plot(x,y1,'*',x,y2,'o',xj,yj1,xj,yj2,xi(2),yi,'rp');

专题五总结

6数值微积分与方程求解

6.1数值微分与数值积分

数值微分：

数值差分：向前差分；向后差分；中心差分；；diff

差商：向前差商；向后差商；中心差商；；

数值积分：牛顿-莱布尼兹公式；梯形法；辛普森法；高斯求积公式；

例1 设f(x)=sinx，在[0，2π]范围内随机采样，计算f’(x)的近似值，并与理论值f’(x)=cosx 进行比较。

1
2
3
4
5
6
\>> x=[0,sort(2*pi*rand(1,5000)),2*pi];
\>> y=sin(x);
\>> f1=diff(y)./diff(x);
\>> f2=cos(x(1:end-1));
\>> plot(x(1:end-1),f1,x(1:end-1),f2);
\>> d=norm(f1-f2)

例2 分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。

1
2
3
4
5
6
7
8
9
\>> format long
\>> f=@(x) 4./(1+x.^2);
\>> [I,n]=quad(f,0,1,1e-8)

\>> [I,n]=quadl(f,0,1,1e-8)

\>> (atan(1)-atan(0))*4

\>> format short

例3 求定积分。
被积函数文件fe.m：

1
2
3
4
5
6
 function f=fe(x)
 f=1./(x.*sqrt(1-log(x).^2));

 

\>> I=integral(@fe,1,exp(1))

例4 求定积分。
被积函数文件fsx.m：

1
2
3
4
5
6
 function f=fsx(x)
 f=sin(1./x)./x.^2;

 

\>> I=quadgk(@fsx,2/pi,+Inf)

例5 设x=1:6，y=[6,8,11,7,5,2],用trapz函数计算定积分。

1
2
3
4
5
6
7
8
 \>> x=1:6;
 \>> y=[6,8,11,7,5,2];
 \>> plot(x,y,'-ko');
 \>> grid on
 \>> axis([1,6,0,11]);
 \>> I1=trapz(x,y)

\>> I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)

例6 分别求二重积分和三重积分。 

1
2
3
4
5
 \>> f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y);
 \>> I1=quad2d(f1,-2,2,-1,1)

\>> f2=@(x,y,z) 4*x.*z.*exp(-z.*z.*y-x.*x);
 \>> I2=integral3(f2,0,pi,0,pi,0,1)

6.2线性方程组求解

直接法：高斯消去法；列主消元法（左除）；矩阵的三角分解法；lu分解

例1 用左除运算符求解下列线性方程组。

1
2
3
\>> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
\>> b=[13,-9,6,0]';
\>> x=A\b

例2 用LU分解求解例1中的线性方程组。

1
2
3
4
\>> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
\>> b=[13,-9,6,0]';
\>> [L,U]=lu(A);
\>> x=U\(L\b)

雅可比迭代法的函数文件jacobi.m：

1
2
3
4
5
6
7
8
9
10
11
12
13
function [y,n]=jacobi(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
  x0=y;
  y=B*x0+f;
  n=n+1;
end

Gauss-Serdel迭代法的函数文件gauseidel.m：

1
2
3
4
5
6
7
8
9
10
11
12
13
function [y,n]=gauseidel(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1); 
U=-triu(A,1);
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=ep
  x0=y;
  y=B*x0+f;
  n=n+1;
end

例3 分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0，迭代精度为10^(-6)。

1
2
3
4
5
 \>> A=[4,-2,-1;-2,4,3;-1,-3,3];
 \>> b=[1,5,0]';
 \>> [x,n]=jacobi(A,b,[0,0,0]',1.0e-6)

\>> [x,n]=gauseidel(A,b,[0,0,0]',1.0e-6)

例4 分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组，看是否收敛。

1
2
3
4
5
 \>> A=[1,2,-2;1,1,1;2,2,1];
 \>> b=[9;7;6];
 \>> [x,n]=jacobi(A,b,[0;0;0],1.0e-6)

\>> [x,n]=gauseidel(A,b,[0;0;0],1.0e-6)

6.3 线性方程组应用举例

1．平面桁架结构受力分析问题

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
alpha=sqrt(2)/2;
A=[0,1,0,0,0,-1,0,0,0,0,0,0,0;
 0,0,1,0,0,0,0,0,0,0,0,0,0;
 alpha,0,0,-1,-alpha,0,0,0,0,0,0,0,0;
 alpha,0,1,0,alpha,0,0,0,0,0,0,0,0;
 0,0,0,1,0,0,0,-1,0,0,0,0,0;
 0,0,0,0,0,0,1,0,0,0,0,0,0;
 0,0,0,0,alpha,1,0,0,-alpha,-1,0,0,0;
 0,0,0,0,alpha,0,1,0,alpha,0,0,0,0;
 0,0,0,0,0,0,0,0,0,1,0,0,-1;
 0,0,0,0,0,0,0,0,0,0,1,0,0;
 0,0,0,0,0,0,0,1,alpha,0,0,-alpha,0;
 0,0,0,0,0,0,0,0,alpha,0,1,alpha,0;
 0,0,0,0,0,0,0,0,0,0,0,alpha,1];
b=[0;10;0;0;0;0;0;15;0;20;0;0;0];
f=A\b;
disp(f')

2．小行星运行轨道计算问题

1
2
3
4
5
6
7
8
xi=[1.02,0.87,0.67,0.44,0.16];
yi=[0.39,0.27,0.18,0.13,0.13];
A=zeros(length(xi));
for i=1:length(xi)
  A(i,:)=[xi(i)*xi(i),2*xi(i)*yi(i),yi(i)*yi(i),2*xi(i),2*yi(i)];
end
b=-ones(length(xi),1);
ai=A\b

或者：

1
2
3
4
5
6
7
 xi=[1.02,0.87,0.67,0.44,0.16]';
 yi=[0.39,0.27,0.18,0.13,0.13]';
 A=[xi.*xi,2*xi.*yi,yi.*yi,2*xi,2*yi];
 b=-ones(length(xi),1);
 ai=A\b
\>> f=@(x,y) 2.4645*x.^2-0.8846*x.*y+6.4917*y.^2-1.3638*x-7.2016*y+1;
 \>> h=ezplot(f,[-0.5,1.2,0,1.2]);

6.4 非线性方程求解与函数极值计算

代法；二分法；割线法

Fzero

Fsolve

求极小值

Fminbnd

Fminsearch

Fminunc

fmincon

例1 求f(x)=0在x0=-5和x0=1作为迭代初值时的根。

1
2
3
4
5
6
7
8
9
10
11
12
13
 \>> f=@(x) x-1./x+5;
 \>> x1=fzero(f,-5)

 \>> x2=fzero(f,1)

 \>> x3=fzero(f,0.1)

\>> f=@(x) x-1./x+5;
 \>> x1=fsolve(f,-5,optimset('Display','off'))

 \>> x2=fsolve(f,1,optimset('Display','off'))

 \>> x3=fsolve(f,0.1,optimset('Display','off'))

f(x)=x^2-1=0的根

1
2
3
4
5
6
7
8
9
10
f=@(x) x.^2-1;
x=[];
x0=-0.25:0.001:0.25;
for x00=x0
  x=[x,fzero(f,x00)];
end
plot(x0,x,'-o')
xlabel('初值');
ylabel('方程的根');
axis([-0.25,0.25,-1,1])

例2 求下列方程组在（1，1，1）附近的解并对结果进行验证。

1
2
3
4
5
6
\>> f=@(x) [sin(x(1))+x(2)+x(3)^2*exp(x(1)),x(1)+x(2)+x(3),x(1)*x(2)*x(3)];
\>> f([1,1,1])

\>> x=fsolve(f,[1,1,1],optimset('Display','off'))

\>> f(x)

例3 求函数f(x)在区间(-10，-1)和(1，10)上的最小值点。

1
2
3
4
\>> f=@(x) x-1./x+5;
\>> [xmin,fmin]=fminbnd(f,-10,-1)

\>> [xmin,fmin]=fminbnd(f,1,10)

例4 求解有约束最优化问题。

1
2
3
4
5
6
7
f=@(x) 0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3;
x0=[0.5;0.5];
A=[-1,-0.5;-0.5,-1];
b=[-0.4;-0.5];
lb=[0;0];
option=optimset('Display','off');
[xmin,fmin]=fmincon(f,x0,A,b,[],[],lb,[],[],option)

例5 要使每周送货车的里程最小，仓库应建在xy平面的什么位置？

1
2
3
4
5
a=[10,30,16.667,0.555,22.2221];
b=[10,50,29,29.888,49.988];
c=[10,18,20,14,25];
f=@(x) sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2));
[xmin,fmin]=fminsearch(f,[15,30])

例6 在例5中，如果由于地域的限制，仓库必须建在曲线y=x^2上，则它应该建在何处？

1
2
3
4
5
6
7
8
9
10
11
12
 非线性约束的函数文件funny.m：
 function [c,h]=funny(x)
 c=[];
 h=x(2)-x(1)^2;

 

\>> a=[10,30,16.667,0.555,22.2221];
 \>> b=[10,50,29,29.888,49.988];
 \>> c=[10,18,20,14,25];
 \>> f=@(x) sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2));
 \>> [xmin,fmin]=fmincon(f,[15,30],[],[],[],[],[],[],'funny')

6.5 常微分方程数值求解

单步法；

多步法

Solver

例1 求解微分方程初值问题，并与精确解进行比较。 

1
2
3
4
f=@(t,y) (y^2-t-2)/4/(t+1);
[t,y]=ode23(f,[0,10],2);
y1=sqrt(t+1)+1;
plot(t,y,'b:',t,y1,'r');

例2 已知一个二阶线性系统的微分方程，取a=2，绘制系统的时间响应曲线和相平面图。

1
2
3
4
f=@(t,x) [-2,0;0,1]*[x(2);x(1)]; 
[t,x]=ode45(f,[0,20],[1,0]); 
subplot(2,2,1);plot(t,x(:,2));
subplot(2,2,2);plot(x(:,2),x(:,1));

例3 假如点燃一个火柴，已知火焰球简化模型，分析λ的大小对方程求解过程的影响。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 lamda=0.01;
 f=@(t,y) y^2-y^3;
 tic;[t,y]=ode45(f,[0,2/lamda],lamda);toc
 disp([‘ode45计算的点数' num2str(length(t))]);

lamda=1e-5;
 f=@(t,y) y^2-y^3;
 tic;[t,y]=ode45(f,[0,2/lamda],lamda);toc
 disp(['ode45计算的点数' num2str(length(t))]);

lamda=1e-5;
 f=@(t,y) y^2-y^3;
 tic;[t,y]=ode15s(f,[0,2/lamda],lamda);toc
 disp(['ode15s计算的点数' num2str(length(t))]);

6.6 常微分方程应用举例

1．Lotka-Volterra模型
第①问：

1
2
3
4
5
6
7
8
rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[300,150])
subplot(1,2,1);plot(t,x(:,1),'-',t,x(:,2),'-*');
legend('x1(t)','x2(t)');
xlabel('时间');ylabel('物种数量');
grid on
subplot(1,2,2);plot(x(:,1),x(:,2))
grid on

第④问：
取λ=0.01， 所以稳定平衡点（1/λ，2/λ）即是（100，200），以此点作为初值时，画出其图像。

1
2
3
4
5
6
7
rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));...
         x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[100,200]);
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');

当将初始值变为（98,195）时，即向下十分接近平衡点，画出其图像。

1
2
3
4
5
6
7
rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));...
         x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[98,195]);
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');

当将初始值变为（70,150）时（向下偏离平衡点比较远时），画出其图像。

1
2
3
4
5
6
7
rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));...
         x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,30],[70,150]);
plot(t,x(:,1),'-o',t,x(:,2),'-*');
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');

当将初始值变为（900,1600）时（向上偏离平衡点十分远时），画出其图像。

1
2
3
4
5
6
7
rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));...
         x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,500],[900,1600]);
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');

  1. Lotka-Volterra改进模型
    第①问：在原模型下，绘制狐狸和兔子数量的时间函数曲线。

    1
    2
    3
    4
    5
    6
    7
    rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));...
             x(2)*(-1+0.01*x(1))];
    plot(t,x(:,1),t,x(:,2));
    legend('x1(t)-兔子','x2(t)-狐狸');
    xlabel('时间');
    ylabel('物种数量');
    title('原模型下，狐狸和兔子数量的函数曲线');

第②问：在改进模型下，狐狸和兔子数量的时间函数曲线。

1
2
3
4
5
6
7
8
rabbitFox=@(t,x) [2*x(1)*(1-x(1)/400-0.005*x(2));...
        x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(t,x(:,1),t,x(:,2));
legend('x1(t)-兔子','x2(t)-狐狸');
xlabel('时间');
ylabel('物种数量');
title('改进模型下，狐狸和兔子数量的函数曲线');

第③问：在原模型下，绘制狐狸数量相对于兔子数量的关系曲线。

1
2
3
4
5
6
7
rabbitFox=@(t,x) [x(1)*(2-0.01*x(2));...
         x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(x(:,1),x(:,2));
xlabel('兔子数量');
ylabel('狐狸数量');
title('原模型下，狐狸数量相对于兔子数量的关系曲线');

第④问：在改进模型下，狐狸数量相对于兔子数量的关系曲线。

1
2
3
4
5
6
7
rabbitFox=@(t,x) [2*x(1)*(1-x(1)/400-0.005*x(2));...
        x(2)*(-1+0.01*x(1))];
[t,x]=ode45(rabbitFox,[0,50],[300,150]);
plot(x(:,1),x(:,2));
xlabel('兔子数量');
ylabel('狐狸数量');
title('改进模型下，狐狸数量相对于兔子数量的关系曲线');

专题六总结

7 matlab 符号计算

7.1 符号对象: symbol

符号计算的结果是一个精确的数学表达式:符号推演，精确解。

数值计算的结果是一个数值：有舍入误差，近似解。

eval(): 将符号表达式转化成数值结果。

syms命令，定义符号变量

assume(condition)

assume(expr,set)

符号对象的运算：

  • 四则运算
  • 关系运算
  • 逻辑运算
  • 因式分解和展开运算：factor(s);expand(s);collect(s);collect(s,v)
  • 其他运算

符号矩阵

梅森素数的验证问题。请验证M19 、 M23 、 M29 、 M31是否为梅森素数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
>> syms p;
\>> m=2^p-1;
\>> p=19;
\>> m19=eval(m)
\>> factor(m19)
\>> p=23;
\>> m23=eval(m)
\>> factor(m23)
\>> p=29;
\>> m29=eval(m)
\>> factor(m29)
\>> p=31;
\>> m31=eval(m)
\>> factor(m31)

例1 求方程ax^2+bx+c=0的根。

1
2
3
4
5
>> syms a b c x;
\>> f=a*x^2+b*x+c
\>> g=coeffs(f,x)
\>> g=g(end:-1:1)
\>> roots(g)

例2 当λ取何值时，以下齐次线性方程组有非零解。

1
2
3
4
syms lamda;
A=[1-lamda,-2,4;2,3-lamda,1;1,1,1-lamda];
D=det(A);
factor(D)

7.2 符号微积分

  1. 极限：limit(f,x,a,’right/left’)

  2. 倒数：diff(f,x,n)

  3. 积分：

    • 不定积分：int(f,x)
    • 定积分：int(f,x,a,b); 无限：inf

例1 求下列极限。

1
2
3
4
5
>> syms a m x n;
\>> f=(x^(1/m)-a^(1/m))/(x-a);
\>> limit(f,x,a)
\>> g=(1+1/n)^n;
\>> limit(g,n,inf)

例2 求下列函数的导数。

1
2
3
4
5
6
7
8
>> syms x y z;
\>> f=sqrt(1+exp(x));
\>> diff(f)  

\>> g=x*exp(y)/y^2;
\>> diff(g,x)
     
\>> diff(g,y)

例3 求下列不定积分。

1
2
3
4
5
>> syms x t;
\>> f=(3-x^2)^3;
\>> int(f)
\>> g=5*x*t/(1+x^2);
\>> int(g,t)

​

例4 求下列定积分。

1
2
3
4
5
6
>> syms x t;
\>> int(abs(1-x),1,2)

\>> int(1/(1+x^2),-inf,inf)

\>> int(4*x/t,t,2,sin(x))

河道水流量的估算问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
xi=0:50:600;
yi=[4.4,4.5,4.6,4.8,4.9,5.1,5.4,5.2,5.5,5.2,4.9,4.8,4.7];
p=polyfit(xi,yi,3);
plot(xi,yi,'o',xi,polyval(p,xi));
syms y x;
y=poly2sym(p,x);
s=int(y,x,0,600);   
v=s*0.6;       
eval(v)

yn=-yi;
p=polyfit(xi,yn,3);
plot(xi,yn,'o',xi,polyval(p,xi));
syms y x yii;
y=poly2sym(p,x);
yii=diff(y,x); 
x=50:60;                  
k=eval(yii);                
all(abs(k)<1/1.5)

7.3 级数

  1. 级数求和
    • sum()
    • symsum(s,v,n,m)

  2. 泰勒级数：将任意函数表示为幂级数

    • taylor(f,v,a,Name,Value)

  3. 复杂函数的计算方法问题

    • 加法= 逻辑电路
    • 减法 = 加法转换
    • 乘除 = 加法 + 移位
    • 复杂函数 = 转换成四则运算，例泰勒级数

例1 求下列级数之和。

1
2
3
4
5
6
>> syms n;
\>> s1=symsum(n^2,1,100)
\>> s2=symsum((-1)^(n-1)/n,1,inf)
\>> s3=symsum((-1)^(n-1)/(2*n-1),n,1,inf)
\>> hypergeom([-1/2, 1], 1/2, -1)   %超几何函数
\>> eval(s3)*4

银行利率的计算问题。

1
2
3
4
5
6
7
8
9
>> syms k r;
\>> p2=symsum(50000*(1+0.045/k)^k,k,2,2);
\>> eval(p2)
\>> p4=symsum(50000*(1+0.045/k)^k,k,4,4);
\>> eval(p4)
\>> p12=symsum(50000*(1+0.045/k)^k,k,12,12);
\>> eval(p12)
\>> limit((1+r/k)^k,k,inf)
\>> 50000*exp(0.045)

例2 求函数f(x)在x=1处的5阶泰勒级数展开式。

1
2
3
4
>> syms x;
\>> f=(1+x+x^2)/(1-x+x^2);
\>> taylor(f,x,1,'Order',6)
\>>expand(ans)

例3 利用泰勒展开式计算三角函数的值。

1
2
3
4
5
>> syms x;
\>> f=taylor(cos(x),x,pi)
\>> x=3;
\>> eval(f)
\>> cos(3)

7.4 符号方程求解

  1. 代数方程：
    • solve(s) ;不一定准确
    • solve(s,v)
    • solve(s1,sn,v1,vn)
  2. 常微分方程 equation;condition;variable
    • D表示倒数: Dy;D2y
    • dsolve(e,c,v)
    • dsolve(e1,e2,c1,cn,v)

例1 解方程ax^2+bx+c=0。 

1
2
3
4
5
6
7
>> syms x y a b c; 
\>> solve(a*x^2+b*x+c==0)
\>> f=a*x^2+b*x+c==0;
\>> solve(f) 
\>> solve(a*x^2+b*x+c)
\>> f=a*x^2+b*x+c
\>> solve(f)

例2 求下列微分方程或方程组的解。 

1
2
3
>> syms x y t;
\>> y=dsolve('Dy-(x^2+y^2)/x^2/2',x)
\>> [x,y]=dsolve('Dx=4*x-2*y','Dy=2*x-y',t)

疾病传染问题。

1
2
3
>> syms a b c y t;
\>> f=dsolve('Dy=a*y*(1-y)-c*y', 'y(0)=b',t)
\>> f=dsolve('Dy=a*y*(1-y)-a*y', 'y(0)=b',t)

专题七知识点总结

8 matlab 图形用户界面设计

8.1 图形窗口与坐标轴

  1. 图形对象：曲线、曲面、坐标轴、图形窗口

    • 句柄

    • 属性

      points 磅

  2. 图形窗口对象

  3. 坐标轴对象

例1 绘制多个图形，并保存图形句柄。

1
2
3
4
5
6
7
8
t=0:pi/10:2*pi;
h1=plot3(t+pi,t-2*pi,sin(t),'r'); 
hold on;
[x,y]=meshgrid(t);
z=sin(x);
h2=mesh(t-2*pi,t+pi,z);
[x3,y3,z3]=cylinder(t);
h3=surf(x3,y3,z3);

例2 分别在两个子图中绘制曲线和曲面。然后设置子图1的背景色为黄色，曲线线条颜色为红色，设置子图2的背景色为青色。

1
2
3
4
5
6
7
8
9
10
subplot(1,2,1)
h1=fplot(@(t)t.*sin(t),@(t)t.*cos(t),[0,6*pi] );
axis equal
subplot(1,2,2)
[x,y,z]=peaks(20);
h2=mesh(x,y,z);
h10=h1.Parent;
h10.Color='y';
h1.Color='r';
h2.Parent.Color='cyan';

例3 建立一个图形窗口。该图形窗口没有菜单条，标题名称为“图形窗口示例”，起始于屏幕左下角、宽度和高度分别为300像素点和150像素点，背景颜色为青色，且当用户从键盘按下任意一个键时，然后在窗口中单击鼠标左键，在鼠标指针所在位置将显示“Hello,World!” 。

1
2
3
4
5
6
7
hf=figure;
hf.Color=[0,1,1];
hf.Position=[1,1,300,150];
hf.Name='图形窗口示例';
hf.NumberTitle='off';
hf.MenuBar='none';
hf.ButtonDownFcn='gtext(''Hello,World!'')';

例4 利用坐标轴对象实现图形窗口的分割。

1
2
3
4
5
6
7
8
ha1=axes('Position',[0.1,0.1,0.7,0.7]);
contour(peaks(20))
ha1.Title=title('等高线');
ha1.YLabel=ylabel('南北向');
ha1.XLabel=xlabel('东西向');
ha2=axes('Position',[0.65,0.7,0.28,0.28]);
surf(peaks(20))
ha2.View=[-30,45];

例5 定义包含4种颜色的ColorOrder属性，绘制6条曲线。

1
2
3
4
5
6
7
x=[0,0];y=[0,1];
ha=axes;
ha.ColorOrder=[0,0,0; 1,0,0; 0,1,0; 0,0,1];
hold on
plot(x,y, x+0.5,y, x+1,y, x+1.5,y, x+2,y, x+2.5,y);
ha.XLim=[-0.2,3];
ha.YLim=[-0.2,1.2];

8.2 曲线与曲面对象

  1. 曲线对象

    • line(x,y,z)

      plot 函数每调用一次，就会刷新坐标轴，清空原有图形再绘制新的曲线

  2. 曲面对象

    • surface(x,y,z,c)

      surf函数每调用一次，就会刷新坐标轴，清空原有图形再绘制新的曲面

  3. 关照处理

    • light（）

  4. 图形对象的反射特性

例1 利用曲线对象绘制五环图案。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
t=-0.1 : 0.1 : 2*pi;
x=cos(t); y=sin(t);
line(x,y,'Color','b')
line(x+1.2,y-1,'Color','y')
line(x+2.4,y,'Color','k')
line(x+3.6,y-1,'Color','g')
line(x+4.8,y,'Color','r')
ha=gca; %获取当前坐标轴句柄
for n=1:size(ha.Children)
  ha.Children(n).LineWidth=5;
end
ha.XLim=[-2,7];
ha.YLim=[-3,2];
axis equal %使圆呈现正圆

例2 利用曲面对象绘制立体圆环

1
2
3
4
5
6
7
8
9
10
11
r=linspace(0,2*pi,60);
[u,v]=meshgrid(r);
x=(8+3*cos(v)).*cos(u);
y=(8+3*cos(v)).*sin(u);
z=3*sin(v);
axes('view',[-37.5,30])
hs=surface(x,y,z);
axis equal;

hs.EdgeColor='none'; %无网格线
hs.FaceColor='interp'; %曲面网格用插值模式填充

例3 绘制光照处理后的圆环面并观察不同光照模式下的效果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
r=linspace(0,2*pi,60); [u,v]=meshgrid(r);
x=(8+3*cos(v)).*cos(u); y=(8+3*cos(v)).*sin(u); z=3*sin(v);
axes('Position',[0.05,0.675,1.0,0.3], 'View',[-37.5,30]);
hs1=surface(x,y,z); axis equal;
hs1.EdgeColor='none'; hs1.FaceColor='interp';
axes('Position',[0.05,0.35,1.0,0.3], 'View',[-37.5,30]);
hs2=surface(x,y,z); axis equal;
hs2.EdgeColor='none'; hs2.FaceColor='interp';
light('Position',[0,0,8])
lighting flat
axes('Position',[0.05,0.025,1.0,0.3], 'View',[-37.5,30]);
hs3=surface(x,y,z); axis equal;
hs3.EdgeColor='none'; hs3.FaceColor='interp';
light('Position',[0,0,8])
lighting phong

例4 绘制具有不同反射特性的圆环面并观察反射特性对图形效果的影响。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
r=linspace(0,2*pi,60); [u,v]=meshgrid(r);
x=(8+3*cos(v)).*cos(u); y=(8+3*cos(v)).*sin(u); z=3*sin(v);
axes('Position',[0.05,0.675,1.0,0.3],'View',[-37.5,30]);
hs1=surface(x,y,z);axis equal;
hs1.EdgeColor='none';hs1.FaceColor='interp';
light('Position',[0,0,8]) ;lighting phong
hs1.SpecularStrength=0.1;
axes('Position', [0.05,0.35,1.0,0.3],'View',[-37.5,30]);
hs2=surface(x,y,z);axis equal;
hs2.EdgeColor='none';hs2.FaceColor='interp';
light('Position',[0,0,8]) ;lighting phong
hs2.SpecularStrength=0.5;
axes('Position', [0.05,0.025,1.0,0.3],'View',[-37.5,30]);
hs3=surface(x,y,z);axis equal;
hs3.EdgeColor='none';hs3.FaceColor='interp';
light('Position',[0,0,8]) ;lighting phong ;
hs3.SpecularStrength=1.0;

8.3 图形用户界面设计方法

Graphical User Interface(GUI)

  1. 控件函数
    • uicontrol()
  2. GUIDE

例1 在图形窗口中建立三个按钮对象，当单击按钮时分别绘制正弦曲线、显示或隐藏坐标轴的网格、清除坐标轴的图形。

1
2
3
4
5
6
7
8
9
10
11
ha= axes('Units','pixels','Position',[40,40,360,360]);  
ptgrid=uicontrol('Style','pushbutton',...
   'String','网格', 'Position', [450,120,48,20],...
  'Callback','grid' );
btncla = uicontrol('Style', 'pushbutton', ...
  'String', '清除','Position', [450,80,48,20],...
  'Callback','cla' );  
btnplot = uicontrol('Style', 'pushbutton', ...
   'String','绘图','Position', [450,160,48,20]); 
%设置“绘图”按钮的Callback属性值是plot_sin函数句柄。
btnplot.Callback=@plot_sin;

新建回调函数文件plot_sin.m

1
2
3
4
function plot_sin(source, callbackdata)
  t=-pi:pi/20:pi;
  plot(t,sin(t));
end

例2 在例1的界面中添加“图形选项”菜单项，其中包括一个二级菜单项“线型”，其下又有3个子菜单项，分别为“实线”、“虚线”、“双划线”。

1
2
3
4
5
6
7
8
hopt=uimenu(gcf,'Label','图形选项','Accelerator','L');
hLStyle=uimenu(hopt,'Label','线型','Tag','LStyle', 'Enable','off');
hL_Solid=uimenu(hLStyle,'Label','实线',...
  'Tag','Solid','Callback', @MLine_Type);
hL_Dotted=uimenu(hLStyle,'Label','虚线',...
  'Tag','Dotted','Callback', @MLine_Type);
hL_Dashed=uimenu(hLStyle,'Label','双划线',...
  'Tag','Dashed','Callback', @MLine_Type);

新建回调函数文件MLine_Type.m

1
2
3
4
5
6
7
8
9
10
function MLine_Type(source,callbackdata)
  hline=findobj('Type','line');
  if strcmp(source.Tag,'Solid' )==1
    hline.LineStyle='-';
  elseif strcmp(source.Tag,'Dotted' )==1
    hline.LineStyle=':'; 
  elseif strcmp(source.Tag,'Dashed' )==1
    hline.LineStyle='--';  
  end
end

修改回调函数文件plot_sin.m

1
2
3
4
5
6
> function plot_sin(source, callbackdata)
>   t=-pi:pi/20:pi;
>   plot(t,sin(t));
>   h1=findobj('Tag','LStyle');
>   h1.Enable='On';  %使得 “线型”菜单项可用
> end

8.4 用户界面设计工具

  1. 回调函数Callback
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%“绘图”按钮的Callback函数
function pushbutton1_Callback(hObject, eventdata, handles)
A=eval(handles.editfz.String);
f=eval(handles.editpl.String)/50;
theta=eval(handles.editxj.String)/180*pi;
x=linspace(0,2*pi,60);
if handles.OpSin.Value==1
  y=A*sin(f*x+theta);
else
  y=A*cos(f*x+theta);
end
plot(x,y);
handles.PStyle.Enable='On';



%“实线”菜单项的回调函数
function Solid_Callback(hObject, eventdata, handles)
hline=findobj('Type','line');
hline.LineStyle='-';
handles.Solid.Checked='On';
handles.Dotted.Checked='Off';
handles.Dashed.Checked='Off';



%“红”菜单项的回调函数
function r_Callback(hObject, eventdata, handles)
hline=findobj('Type','line');
hline.Color='r';
handles.r.Checked='On';
handles.g.Checked='Off';
handles.b.Checked='Off';

8.5 APP 设计工具

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
%绘图按钮的回调函数
function Draw_Plot(app)
f=eval(app.DiscreteKnob.Value); %获取离散旋钮对象的Value属性值
theta=app.Knob.Value/180*pi; %获取旋钮对象的Value属性值
x=linspace(0,2*pi,60);
if app.RadioButton.Value==1 
y=sin(f*x+theta);
else
y=cos(f*x+theta);
end
% plot函数的第1个参数必须是坐标轴对象句柄
%plot函数默认的线条宽度为1，App程序的线条宽度必须小于0.4
plot(app.UIAxes,x,y,'LineWidth',0.2);
end



%清空按钮的回调函数
function Clear_Axes(app)
cla(app.UIAxes)
end


8.6 图形用户界面应用举例

例1 利用GUIDE设计工具设计如图所示的用户界面。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
%图形窗口的打开事件响应
handles.peaks=peaks(34);
handles.membrane=membrane;
[x,y]=meshgrid(-8:0.3:8);
r=sqrt(x.^2+y.^2);
sinc=sin(r)./(r+eps);
handles.sinc=sinc;
handles.current_data=handles.sinc;
colormap(spring);

%Mseh_Callback函数
    mesh(handles.current_data)

%Surf_Callback函数
    surf(handles.current_data)

%Contour3_Callback函数
    contour3(handles.current_data)

%切换按钮的Callback函数
if hObject.Value==1
  grid on
  hObject.String='隐藏网格';
else 
  grid off
  hObject.String='显示网格';
end

%列表框的Callback函数
str=hObject.String;
val=hObject.Value;
switch strtrim(str{val})
case 'Peaks' 
  handles.current_data=handles.peaks;
case 'Membrane' 
  handles.current_data=handles.membrane;
case 'Sinc'
  handles.current_data=handles.sinc;
end
guidata(hObject,handles)

%色图弹出式菜单的Callback函数
str=hObject.String;
cm=hObject.Value;
colormap( eval(str{cm}) );

%视点设置按钮的Callback函数
el=eval(handles.edit_el.String);
az=eval(handles.edit_az.String);
view(az,el)

%着色方式按钮组ChooseShading的SelectionChanged函数
switch eventdata.NewValue.Tag
  case 'rb_flat'
    shading flat;
  case 'rb_interp'
    shading interp;
  case 'rb_faceted'
    shading faceted;
end



例2 生成一个用于观察周期信号波形叠加效果的程序模块。

1
2
3
4
5
6
7
8
9
10
%“开始”按钮的回调函数
amr=eval(app.ARadio.Value);
theta=app.PhDiff.Value/180*pi;
x=linspace(0,10*pi,300);
y=sin(x)+sin(3*x+theta)/amr;
if strcmp(app.Noise.Value,'On')
        y=awgn(y,30);
end
app.WaveA.Value=max(y);
plot(app.UIAxes,x,y,'LineWidth',0.1);

专题八总结

Simulink 系统仿真

启动

  1. “主页”选项卡：新建
  2. “主页”选项卡：simulink
  3. 命令行窗口，simulink

模型

模块

模型存盘

模块参数

步骤

  1. 建立系统仿真模型
  2. 设置仿真参数
  3. 启动仿真并分析仿真结果

9.2 子系统的创建与封装

创建

  1. 新建
  2. 模块转换

封装

条件执行

9.3 s函数的设计与应用

S函数： 系统函数，system function

采用S函数实现y=kx+b。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
%①主函数
function [sys,x0,str,ts]=timekb(t,x,u,flag,k,b)
switch flag
  case 0
    [sys,x0,str,ts]=mdlInitializeSizes; %初始化
  case 3
    sys=mdlOutputs(t,x,u,k,b);     %计算输出量
  case {1,2,4,9}
    sys=[];
  otherwise                %出错处理
    error(num2str(flag))
end


%②初始化子函数
function [sys,x0,str,ts]=mdlInitializeSizes()
sizes=simsizes;
sizes.NumContStates=0;  %无连续状态
sizes.NumDiscStates=0;   %无离散状态
sizes.NumOutputs=1;       %有一个输出量
sizes.NumInputs=1;     %有一个输入信号
sizes.DirFeedthrough=1;  %输出量中含有输入量
sizes.NumSampleTimes=1;  %单个采样周期
sys=simsizes(sizes);
%给其他返回参数赋值
x0=[];          %设置初始状态为零状态
str=[];          %将str变量设置为空字符串
ts=[-1,0];         %假定继承输入信号的采样周期


%③输出子函数
function sys=mdlOutputs(t,x,u,k,b)
sys=k*u+b;

9.4 simulink仿真应用举例

仿真求解：

（1）求最大安全体重

1
2
3
4
5
6
7
8
9
10
11
for m=100:-0.5:20
  [t,x,y_w]=sim('bengji',0:0.01:100);
  if min(y_w)>1.5
    break;
  end
end
disp(['最大安全体重是',num2str(m)])
dis=min(y_w);
disp(['最小的安全距离是',num2str(dis)])
plot(t,y_w)
grid on

（2）求最小弹性系数

1
2
3
4
5
6
7
8
9
10
11
12
m=65; 
for k=10:0.1:50
  [t,x,y_w]=sim('bengji',0:0.01:100);
  if min(y_w)>1.5
    break;
  end 
end
disp(['最小弹性系数k是',num2str(k)])
dis=min(y_w);
disp(['最小的安全距离是',num2str(dis)])
plot(t,y_w)
grid on

专题九总结

10 外部程序接口

10.1 在Excel中使用matlab

Spreadsheet Link

10.2 matlab 文件操作

  1. fid = fopen(filename,permission)
    • fid, 文件识别号，不成功时为-1
    • permission，’r’，只读；’w’,写；’a’,追加；‘r+’,读写。
  2. status=fclose(fid)
    • fid,文件标识号，为all时，关闭已所有打开的文件
    • status，0，关闭成功；-1，关闭不成功
  3. 读写文本文件
    • [A,count]=fscanf(fid,fmt,size)
    • count=fprintf(fid,fmt,A)
  4. 读写二进制文件
    • [A,count]=fread(fid,size,precision,skip)
    • count=fwrite(fid,A,precious)
  5. 数据文件定位：文件位置指针
    • fseek(fid,offset,origin)
    • ftell(fid)
    • status=feof(fid)

例1 文件“观测记录.txt”的内容如图所示，读取文件前10行的数据。

1
2
3
4
5
6
7
8
9
10
11
12
fid=fopen('观测记录.txt','r');
title=fscanf(fid,'%s',6);
qxsj=[];
for i=1:10
  qxsj{i,1}=fscanf(fid,'%s',1);
  qxsj{i,2}=fscanf(fid,'%s',1);
  qxsj{i,3}=fscanf(fid,'%f',1);
  qxsj{i,4}=fscanf(fid,'%f',1);
  qxsj{i,5}=fscanf(fid,'%f',1);
  qxsj{i,6}=fscanf(fid,'%s',1);
end
fclose(fid);

例2 计算y=exp(x)sinx，其中x∈[0,2π]。将x、y写入二进制文件“模拟数据.dat”。

1
2
3
4
5
fid=fopen('模拟数据.dat','w');
x=linspace(0,2*pi,100);
y=exp(x).*sin(x);
count=fwrite(fid, [x; y], 'double');
fclose(fid);

例3 读取例2生成的文件中的后40组数据，并绘制图形。

1
2
3
4
5
6
7
8
9
fid=fopen('模拟数据.dat','r');
status=fseek(fid, -40*2*8, 'eof');
x=[]; y=[];
while ~feof(fid)
  x=[x; fread(fid,1,'double')];
  y=[y; fread(fid,1,'double')];
end
plot(x,y);
fclose(fid);

10.3 在其他语言程序中读写matlab的数据文件

接口函数

mat文件：matlab存储数据的标准格式；变量的名称、类型和值均保存。

  • save
  • load
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
#include "mat.h"
\#include <iostream>
using namespace std;

int main()
{
 MATFile *pmat; /* 定义MAT文件指针*/
 mxArray *pa1, *pa2, *pa3;
  //定义要写入的数据
 double data[9] = {1.1, 2.2, 3.3, 4.4, 5.5, 6.6, 7.7, 8.8, 9.9};
 const char *file = "mattest.mat";//定义要操作的文件
 int status; //定义存储函数返回值的变量
 /* 打开一个MAT文件，如果不存在则创建一个MAT文件，如果打开失败，则返回 */
 cout << "生成文件 : " << file << endl;
 pmat = matOpen(file, "w");
 if (pmat == NULL) {
 cout << "不能创建文件 : " << file << endl;
 cout << "(请确认是否有权限访问指定文件夹?)\n";
 return(0);
 }
 /* 创建3个mxArray对象，其中pa1存储一个实数，pa2为3×3的双精度实型矩阵，*/
 /* pa3存储字符串，如果创建失败则返回 */
 pa1 = mxCreateDoubleScalar(1.234);
 if (pa1 == NULL) {
 cout << "不能创建变量。\n";
 return(0);
 }
 pa2 = mxCreateDoubleMatrix(3, 3, mxREAL);
 if (pa2 == NULL) {
 cout << "不能创建矩阵。\n";
 return(0);
 }
 memcpy((void *)(mxGetPr(pa2)), (void *)data, sizeof(data));
 pa3 = mxCreateString("MAT文件示例");
 if (pa3 == NULL) {
 cout << "不能创建字符串。\n";
 return(0);
 }
 /* 向MAT文件中写数据，失败则返回 */
 status = matPutVariable(pmat, "LocalDouble", pa1);
 if (status != 0) {
 cout << "写入局部变量时发生错误。\n " ;
 return(0);
 }
 status = matPutVariableAsGlobal(pmat, "GlobalDouble", pa2);
 if (status != 0) {
 cout << "写入全局变量时发生错误。\n";
 return(0);
 }
 status = matPutVariable(pmat, "LocalString", pa3);
 if (status != 0) {
 cout << ": 写入String类型数据时发生错误。\n ";
 return(0);
 }
 /* 清除矩阵 */
 mxDestroyArray(pa1);
 mxDestroyArray(pa2);
 mxDestroyArray(pa3);
 /* 关闭MAT文件 */
 if (matClose(pmat) != 0) {
 cout << "关闭文件时发生错误 " << endl;
 return(0);
 }
 cout << "文件创建成功！\n";
 return(1);
}

10.4 在matlab中调用其他语言编写的程序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
#include "mex.h"
bool isCoprime(double *inx,double *outy) 
{ 
    int x=*inx,y=*outy,tmp;
    if (x<y) 
    {tmp=x; x=y; y=tmp;}
    tmp=x%y; 
    while (tmp) 
    { 
       x=y; y=tmp; tmp=x%y; 
    } 
    if (y==1)     //最大公约数为1,所以互质 
      return true; 
    else       //最大公约数大于1，所以不互质 
      return false; 
} 

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] )
{
 double *result;
 int m,n,i;
 //检查调用时输入、输出是否符合要求
 if(nrhs!=2) {
   mexErrMsgTxt("输入参数应有两个！"); return; }
  if(nlhs!=1) {
   mexErrMsgTxt("输出参数应只有1个！"); return; }
  for(i=0;i<2;i++){
   m = int(mxGetM(prhs[i]));
   n = int(mxGetN(prhs[i]));
   if( !mxIsDouble(prhs[i]) || m!=1 || n!=1) {
      mexErrMsgTxt("输入参数必须是单个的数！"); return; }
   }
 //在MATLAB工作区建立一个矩阵，用于存放计算结果，
 //矩阵的大小与输入实参的大小相同。 
 plhs[0] = mxCreateDoubleMatrix(m,n, mxREAL);
 result = mxGetPr(plhs[0]);
 if (isCoprime(mxGetPr(prhs[0]),mxGetPr(prhs[1])))
   *result=1;
 else
   *result=0;
}

10.5 在其他语言程序中调用matlab函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <engine.h>
\#include <iostream>
using namespace std;

int main()
{
  //定义一个engine类型的指针ep，用于操作引擎对象。
  Engine *ep; 
  //定义mxArray 类型的指针，用于指向所调用MATLAB函数的输入对象。
  mxArray *T = NULL, *R = NULL;
  //定义存储数据的变量，并按给定表达式生成数据。
  double t[180],r[180];
  double a,b;
  a=2;b=3;
  for (int i=0;i<180;i++)  {
    t[i]=i*0.1;
    r[i]=a+b*t[i]; }
  //启动MATLAB计算引擎
  //Windows系统中，engOpen函数的参数为空字符串，指定在本机启动计算引擎。 
  if (!(ep = engOpen("")))  {
    cout<< "\n不能启动MATLAB引擎\n"; 
    return 0; 
    }
  //建立MATLAB变量，调用MATLAB函数
  T = mxCreateDoubleMatrix(1, 180, mxREAL);
  memcpy((void *)mxGetPr(T), (void *)t, sizeof(t));
  engPutVariable(ep, "T", T); 
  R = mxCreateDoubleMatrix(1, 180, mxREAL);
  memcpy((void *)mxGetPr(R), (void *)r, sizeof(r));
  engPutVariable(ep, "R", R); 
  engEvalString(ep, "polar(T,R);");
  engEvalString(ep, "title('阿基米德螺线');");  
  //释放内存空间，关闭计算引擎
  system("pause");
  mxDestroyArray(T);
  mxDestroyArray(R);
  engClose(ep);
  return 1; 
}

专题十总结

其他学习资源

https://matlabacademy.mathworks.com/cn

https://ww2.mathworks.cn/videos/getting-started-with-matlab-68985.html

发表于 更新于
本文字数： 701 阅读时长 ≈ 1 分钟

技巧积累

[toc]

20200623 caj转pdf

  1. 虚拟打印
  • cajviewer 打印输出
  • 福昕浏览器 打印机
  • 微软打印机
  • 只能单篇
  1. 在线转换工具
  • speedPDF
  • 迅捷PDF转化器
  • 免费次数和大小有限
  1. 知网直接下载PDF
  1. 转换软件
  1. 直接修改后缀名？

20200624 github空间大小

  1. GitHub 有空间限制吗？

没有限制。推荐1G以内。达到du1G以后会受到GITHUB的通知邮件，上传超过50M的单个文件会zhiwarning，无法上传超过100M的单个文件。dao目前大文件会提供一个1G的免费GIT-LFS空间。

  1. 如何查看GitHub空间使用情况？

  2. github更改密码之后，Gitnote还能自动同步吗？

不能，对Gitnote而言，更改GitHub密码要慎重。

解决办法？重新将仓库复制到本地。

对git而言则无影响。

Gitnote仅在克隆仓库时需要输入GitHub账户和密码，这意味着，GitHub其实可以在本地同时为多个GitHub提供同步服务。

图床呢？还能继续放到GitHub吗？

title

验证发现，是可以的。

20200626 manim—用于数学科普，可视化展示。

基本原理：绘图得到图像帧，压缩成视频。

使用思路：场景类（Scene），继承

代码，不同层次上的封装。