先发代码为敬
1、牛顿下山法
#下山法
from sympy import*
x=Symbol('x')
c=Symbol('c')
print('请输入形如a*exp(1)**(b*x)+c+d*x+e*x**2+f*x**3+g*x**4=0的各参数值')
a,b,c,d,e,f,g=eval(input('(a,b,c,d,e,f,g)='))
y=a*x*exp(1)**(b*x)+c+d*x+e*x**2+f*x**3+g*x**4'''
dy=diff(y,x)
z=x-y/c
while True:
ed=str(input('请回车继续，end结束'))
if ed=='end':
break
i=N((input('请输入初始值x0=')),50)
c=dy.subs('x',i)
i1=z.subs([('x',i),('c',c)])
r=int(input('请输入有效数字位数大于多少位,默认大于5位:')or(5))
err=N(10**(-r),50)
j=0
print('x'+str(j+1)+'=',N(i1,50))
while abs(N((i1-i),50))>err:
j=j+1
i1=N(i,50)
i=N(z.subs([('x',N(i,50)),('c',c)]),50)
print('x'+str(j+1)+'=',N(i,50))
print('end')
上周已经使用过sympy解牛顿法。
这里总结一下sympy里的subs方法的用法：
subs()可以替换表达式中的符号，它有如下3种调用方式：
-
expression.subs(x,y)：将算式中的x替换为y； -
expression.subs({x:y,u:v})：使用字典进行多次替换； -
expression.subs([(x,y),(u,v)])：使用列表进行多次替换。
多次替换是顺序执行的，因此expression.subs([(x,y),(u,v)])并不能对符号x和y进行交换。
ps：此次便用了多次替换。
2、循环嵌套及取余取整
#取余。取一个任意小于1美元的金额，然后计算可以换成最少多少枚硬币。硬币有1美分、5美分、10美分、25美分4种。1美元等于100美分。举例来说，0.76美元计算结果应该是3枚25美分，1枚1美分。类似于76枚1美分，2枚25美分+2枚10美分+1枚5美分+1枚1美分这样的结果都是不符合要求的。
x=int(input('请输入小于100的整数'))#x小于100
for a in range(0,x):
for b in range(0,x):
for c in range(0,x):
for d in range(0,x):
if d==(x-a*25-b*10-c*5):
y=[a,b,c,d]
i=0
for i in range(0,4):
if y[i]!=0:
#至少两项非零
if ((y[i-1]==0)&(y[i-2]==0) and (y[i-3]==0))!=True:
#至少有一项为0
if ((y[i-1]!=0)&(y[i-2]!=0) and (y[i-3]!=0))!=True:
print(str(y[0])+'枚25美分',str(y[1])+'枚10美分',str(y[2])+'枚5美分',str(y[3])+'枚1美分')
- 通常用for循环遍历
- while循环判断为真循环继续
- math.ceil向上取整，标准库里//取整，%取余
- 位运算符与布尔运算符
下表中变量 a 为 60，b 为 13，二进制格式如下：
a = 0011 1100
b = 0000 1101
a&b = 0000 1100
a|b = 0011 1101
a^b = 0011 0001
~a = 1100 0011
位运算符
以下假设变量 a 为 10, b为 20:
布尔运算符
3、简单回归计算（相关性系数和纳什效率系数）
import time
import numpy as np
time_start=time.time()
a=np.loadtxt('1.csv',delimiter=',')
i=np.shape(a)[1]
xx=np.mean(a[0])
yy=np.mean(a[1])
C=np.var(a[0])*i
D=np.var(a[1])*i
B=E=F=G=0
for j in range(0,i):
x=a[0][j]
y=a[1][j]
b=(x-y)**2
e=(x-xx)
f=(y-yy)
g=e*f
B=B+b
E=E+e
F=F+f
G=G+g
NS=1-B/C
R2=(G**2)/(C*D)
time_end=time.time()
print(i)
print('纳什效率系数NS等于'+str(NS))
print('相关性系数R2等于'+str(R2))
print('总共耗时'+str(time_end-time_start)+'秒')
此处是numpy的简单应用
国际按例，以下方式导入NumPy库：
import numpy as np
NumPy使用ndarray对象表示数组，它是整个库的核心对象。
如下创建多维数组：
a=np.array([[1,2],[3,4]])
数组的形状可以由shape属性获得，
a.shape
或者方法
np.shape(a)
此处输出为(2,2)
可以用loadtxt()从数据文件载入数据，官方说明如下
numpy.loadtxt(fname, dtype=<type 'float'>, comments='#', delimiter=None, converters=None, skiprows=0, usecols=None, unpack=False, ndmin=0)
- delimiter : str, optional
The string used to separate values. By default, this is any whitespace.
4、简单绘图
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
x = np.linspace(-10, 10, 10000)
y=x*np.exp(1)**x-np.exp(2)**x
z=x**3-x*6
plt.plot(x,y,label="$y$",color="red",linewidth=1)
plt.plot(x,z,'b--',label='$z$',color='blue',linewidth=1)
plt.xlabel("x")
plt.ylabel("y,z")
plt.title(u'函数图')
plt.ylim(-10,10)
plt.legend()
plt.show()
引入简单绘图的包
import matplotlib.pyplot as plt
最后用plt.show()显示图像
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
基本画图方法：plt.plot(x,y)，plt.xlabel(‘x’)，plt.ylabel(‘y’)，plt.title(‘…’)
子图：plt.subplot(abc)，其中abc分别表示子图行数、列数、序号
创建绘图组件的顶层容器：fig = plt.figure()
添加子图：ax = fig.add_subplot(abc)
设置横轴上的主定位器：ax.xaxis.set_major_locator(…)
绘制方图：plt.hist(a,b)，a为长方形的左横坐标值，b为柱高
绘制散点图：plt.scatter(x,y,c = ‘..’,s = ..)，c表示颜色，s表示大小
添加网格线：plt.grid(True)
添加注释：如ax.annotate('x', xy=xpoint, textcoords='offsetpoints',xytext=(-50, 30), arrowprops=dict(arrowstyle="->"))
增加图例：如plt.legend(loc='best', fancybox=True) 对坐标取对数：横坐标plt.semilogx()，纵坐标plt.semilogy()，横纵同时plt.loglog()
目前在读，白话深度学习与TensorFlow，正在考虑将研究方向BMPs与相关算法结合。
以上，下周见。
网友评论