数模竞赛-Python基础

Python入门

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
from numpy.random import randint
import numpy
a=randint(10,20,16)
ma=max(a)
ind2=numpy.where(a==ma)
print(ind2[0])

import string,random,collections
x=string.ascii_letters+string.digits
y=''.join([random.choice(x)for i in range(1000)])
count=collections.Counter(y)
for k,v in sorted(count.items()):
    print(k,':',v)

import random
x=random.randint(1e5,1e8)
y=list(map(int,str(x)))
z=list(map(lambda x,y:x%2==1 and y%2==0,[1,3,2,4,1],[3,2,1,2]))

a=filter(lambda x:x>10,[1,11,2,45,7,6,13])
b=filter(lambda x:x.isalnum(),['abc','xy12','***'])

基础操作

NumPy

arange()生成$[\mathrm{start},\mathrm{end})$上步长间隔为$\mathrm{step}$的等差数组。linspace()生成$[\mathrm{start},\mathrm{end}]$上间隔相等的$\mathrm{num}$个数据的等差数组。logspace()生成区间$\left[\mathrm{base}^{\mathrm{start}},\mathrm{base}^{\mathrm{end}}\right]$上$\mathrm{num}$个数据的等比数组。

1
2
3
4
5
6
7
8
9
10
11
import numpy
a1=numpy.array([1,2,3,4])
a2=a1.astype(float)
a3=numpy.array([1,2,3,4],dtype=float)
b=numpy.array([[1,2,3],[4,5,6]])
c=numpy.arange(1,5) #[1,2,3,4]
    #arange(start=None,end=None,step=None,dtype=None)
d=numpy.linspace(1,4,4) #[1,2,3,4]
    #linspace(start,stop,num=50,endpoint=True)
e=numpy.logspace(1,3,3,base=2) #[2,4,8]
    #logspace(start,end,num=50,endpoint=True,base=10.0)

特殊函数：

1
2
3
4
5
6
7
8
9
import numpy
a=numpy.ones(4,dtype=int) #[1,1,1,1]
b=numpy.ones((4,),dtype=int) #[1,1,1,1]
c=numpy.ones((4,1)) #4行1列
d=numpy.zeros(4) #[0,0,0,0]
e=numpy.empty(3)
f=numpy.eye(3) #3阶单位阵
g=numpy.eye(3,k=1) #第k对角线为1 其他元素为0
h=numpy.zeros_like(a) #与a同维全0数组

元素索引：

1
2
3
4
5
6
7
import numpy
a=numpy.arange(16).reshape(4,4) #生成4行4列
b=a[1][2] #6
c=a[1,2] #6
d=a[1:2,2:3] #[[6]]
x=numpy.array([0,1,2,1])
print(a[x==1]) #a的2、4行元素

vstack()上下合并矩阵，hstack()左右合并：

1
2
3
4
5
6
import numpy
a=numpy.arange(16).reshape(4,4)
b=numpy.floor(5*numpy.random.random((2,4)))
c=numpy.ceil(6*numpy.random.random((4,2)))
d=numpy.vstack([a,b])
e=numpy.hstack([a,c])

vsplit与hsplit分割矩阵：

1
2
3
4
5
import numpy
a=numpy.arange(16).reshape(4,4)
b=numpy.vsplit(a,2) #平均分为2行
print(b[0],b[1])
c=numpy.hsplit(a,4) #平均分为4列

矩阵求和：

1
2
3
4
5
6
7
import numpy
a=numpy.array([[0,3,4],[1,6,4]])
b=a.sum() #所有元素的和
c1=sum(a) #逐列元素的和
c2=numpy.sum(a,axis=0) #同上
c3=numpy.sum(a,axis=0,keepdims=True) #同上
print(c2.shape,c2.shape) #(3,)和(1,3)

+、-、*、/、**都是对应逐个元素运算，这里略。

矩阵乘法：

1
2
3
4
import numpy
a=numpy.ones(4)
b=numpy.arange(2,10,2)
c=a@b #a行向量 b列向量

计算范数：

1
2
3
4
5
6
7
8
9
import numpy
a=numpy.array([[0,3,4],[1,6,4]])
b=numpy.linalg.norm(a,axis=1) #行向量2范数
    #norm(x,ord=None,axis=None,keepdims=False)
c=numpy.linalg.norm(a,axis=0) #列向量2范数
d=numpy.linalg.norm(a) #矩阵2范数
print(numpy.round(b,4))
print(numpy.round(c,4))
print(round(d,4))

求线性方程组唯一解：

$$
\begin{cases}
\displaylines{
3x+y=9,\\
x+2y=8
}
\end{cases}
$$

1
2
3
4
5
import numpy
a=numpy.array([[3,1],[1,2]])
b=numpy.array([9,8])
x1=numpy.linalg.inv(a)@b #第一种方法
x2=numpy.linalg.solve(a,b) #第二种方法

求超定线性方程组最小儿乘解：

$$
\begin{cases}
\displaylines{
3x+y=9,\\
x+2y=8,\\
x+y=6
}
\end{cases}
$$

1
2
3
4
5
import numpy
a=numpy.array([[3,1],[1,2],[1,1]])
b=numpy.array([9,8,6])
x=numpy.linalg.pinv(a)@b
print(numpy.round(x,4))

矩阵特征值和特征向量：

$$
\left[
\begin{matrix}
\displaylines{
0&0&0&1\\
0&0&1&0\\
0&1&0&0\\
1&0&0&0
}
\end{matrix}
\right]
$$

1
2
3
4
import numpy
a=numpy.eye(4)
b=numpy.rot90(a)
c,d=numpy.linalg.eig(b) #特征值 特征向量

SciPy

求方程$x^{980}-5.01x^{979}+7.398x^{978}-3.388x^{977}-x^3+5.01x^2-7.398x+3.388=0$在给定初值$1.5$附近的一个实根。

1
2
3
4
from scipy.optimize import fsolve,root
fx=lambda x:x**980-5.01*x**979+7.398*x**978-3.388*x**977-x**3+5.01*x**2-7.398*x+3.388
x1=fsolve(fx,1.5,maxfev=4000) #法一
x2=root(fx,1.5) #法二

求方程组的一组数值解：

$$
\begin{cases}
\displaylines{
x_1^2+x_2^2=1,\\
x_1=x_2
}
\end{cases}
$$

1
2
3
4
from scipy.optimize import fsolve,root
fx=lambda x:[x[0]**2+x[1]**2-1,x[0]-x[1]]
s1=fsolve(fx,[1,1])
s2=root(fx,[1,1])

分别计算$\begin{cases}\displaylines{a=2,\\b=1}\end{cases}$和$\begin{cases}\displaylines{a=2,\\b=10}\end{cases}$时，$\displaystyle\int_0^1\left(ax^2+bx\right)\mathrm dx$的值。

1
2
3
4
5
from scipy.integrate import quad
def fun46(x,a,b):
    return a*x**2+b*x
I1=quad(fun46,0,1,args=(2,1))
I2=quad(fun46,0,1,args=(2,10))

最小二乘解略。

求矩阵最大模特征值及对应的特征向量：

$$
\left[
\begin{matrix}
\displaylines{
1&2&3\\
2&1&3\\
3&3&6
}
\end{matrix}
\right]
$$

1
2
3
4
5
from scipy.sparse.linalg import eigs
import numpy
a=numpy.array([[1,2,3],[2,1,3],[3,3,6]],dtype=float)
b,c=numpy.linalg.eig(a)
d,e=eigs(a,1) #最大模特征值 对应特征向量

SymPy

求代数方程的解：

$$
ax^2+bx+c=0
$$

1
2
3
4
import sympy
a,b,c,x=sympy.symbols('a b c x')
x0=sympy.solve(a*x**2+b*x+c,x)
print(x0)

求方程组的符号解：

$$
\begin{cases}
\displaylines{
x_1^2+x_2^2=1,\\
x_1=x_2
}
\end{cases}
$$

1
2
3
4
5
6
7
8
9
10
11
import sympy

#法一
sympy.var('x1 x2')
s=sympy.solve([x1**2+x2**2-1,x1-x2],[x1,x2])

#法二
x=sympy.var('x:2')
s=sympy.solve([x[0]**2+x[1]**2-1,x[0]-x[1]],x)

print(s)

求矩阵特征值和特征向量的符号解：

$$
\left[
\begin{matrix}
\displaylines{
0&0&0&1\\
0&0&1&0\\
0&1&0&0\\
1&0&0&0
}
\end{matrix}
\right]
$$

1
2
3
4
import sympy
a=sympy.Matrix([[0,0,0,1],[0,0,1,0],[0,1,0,0],[1,0,0,0]])
print(a.eigenvals()) #特征值
print(a.eigenvects()) #特征向量

Matplotlib

常用属性：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
b 蓝色
c 青色
g 绿色
k 黑色
m 品红
r 红色
w 白色
y 黄色
- 实线
-- 虚线
-. 点划线
: 点线
. 点
o 圆圈
* 星型
x 十字架
s 正方形
p 五角星
D/d 钻石/小钻石
h 六角形
+ 加号
| 竖直线
V^<> 下三角 上三角 左三角 右三角
1234 Tripod向下、上、左、右

折线图：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import pandas,pylab
pylab.rc('font',family='SimHei') #中文标签
pylab.rc('font',size=16) #字体大小
pylab.rc('axes',unicode_minus=False) #正常显示负号
a=pandas.read_excel("data.xlsx",header=None)
b=a.values
x=b[0]
y=b[1:]
pylab.plot(x,y[0],'-*b',label='钻石')
pylab.plot(x,y[1],'--dr',label='铂金')
pylab.xlabel('月份')
pylab.ylabel('每月销量')
pylab.legend(loc='upper left')
pylab.grid()
pylab.show()

柱状图：

1
2
3
4
5
6
7
8
9
10
import pandas,pylab
pylab.rc('font',family='SimHei')
pylab.rc('font',size=16)
a=pandas.read_excel("data.xlsx",header=None)
b=a.T
b.plot(kind='bar')
pylab.legend(['钻石','铂金'])
pylab.xticks(range(5),b[0],rotation=0)
pylab.ylabel('数量')
pylab.show()

子图+$\LaTeX$渲染：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#需要apt安装 texlive-latex-base texlive-latex-extra cm-super dvipng
import pylab,numpy
pylab.rc('text',usetex=True)

pylab.subplot(221) #整个图分成2行2列 序号为1 左上角
y1=numpy.random.randint(2,5,6)
y1=y1/sum(y1)
s=["Apple","grape","peach","pear","banana","pineapple"]
pylab.barh(s,y1) #横向柱状图

pylab.subplot(222) #整个图分为2行2列 序号为2 右上角
pylab.pie(y1,labels=s) #饼图

pylab.subplot(212) #整个图分为2行1列 序号为2 正下方
x2=numpy.linspace(0.01,10,100)
y2=numpy.sin(10*x2)/x2
pylab.plot(x2,y2) #折/线图
pylab.xlabel('$x$')
pylab.ylabel('$\\mathrm{sin}(10x)/x$')
pylab.show()

三维曲线,绘制$\begin{cases}\displaylines{x=s^2\sin s,\\y=s^2\cos s,\\ z=s,}\end{cases}s\in[-50,50]$：

1
2
3
4
5
6
7
import pylab,numpy
ax=pylab.axes(projection='3d') #设置为三维图形模式
z=numpy.linspace(-50,50,1000)
x=z**2*numpy.sin(z)
y=z**2*numpy.cos(z)
ax.plot(x,y,z,'k')
pylab.show()

绘制三维表面图$z=50\sin(x+y)$：

1
2
3
4
5
6
7
import pylab,numpy
x=numpy.linspace(-4,4,100)
x,y=numpy.meshgrid(x,x) #x,y轴分别平铺
z=50*numpy.sin(x+y)
ax=pylab.axes(projection='3d')
ax.plot_surface(x,y,z,color='y')
pylab.show()

同样三维表面图，带冷暖图例，绘制$z=\sin(\sqrt{x^2+y^2})$：

1
2
3
4
5
6
7
8
9
import pylab,numpy
ax=pylab.axes(projection='3d')
x=numpy.arange(-6,6,0.25)
y=numpy.arange(-6,6,0.25)
x,y=numpy.meshgrid(x,y)
z=numpy.sin(numpy.sqrt(x**2+y**2))
surf=ax.plot_surface(x,y,z,cmap='coolwarm')
pylab.colorbar(surf)
pylab.show()