python科学计算的⼏个例⼦python解常微分⽅程
python解常微分⽅程的步骤如下：
1. 将计算区间分为n个⼩段，在每⼀⼩段上将求解的曲线作为直线处理；
2. 将⼀个n阶常微分⽅程转换成[y_n,y_n-1,…,y_i,…,y_0]向量的线性⽅程组，其中y_i表⽰y的i阶导数；
3. 确定初始状态
4. 迭代求解每⼀个点的y值（欧拉法），最后通过matplotlib做出曲线图。
以下⾯的三阶常系数微分⽅程为例：
1. 假设求解区间为0~4，划分为400个⼩区间
2. 转换成对应的线性⽅程组，如下图
3. 初始状态设为y2_0 = 2， y1_0 = 1， y0_0 = 0
4. 求解及出图，见代码
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt
y2_0 =2
y1_0 =1
y0_0 =0
y0 =[y0_0, y1_0,y2_0]
# 函数func，y是指上⾯提到的有y的各阶导数组成的向量，t是指⾃变量
def func(y, t):
T=[[1,1,1],[1,0,0],[0,1,0]]
return T@y
x = np.arange(0,4.0,0.1)
t=x
y = odeint(func, y0, t)
print(x)
print(y)
plt.plot(x,y[:,0])
plt.show()
结果如下：
python求解点绕定轴转动后的坐标
⽤python计算定点绕定轴转动⼀定⾓度的坐标，输出结果如下：
import numpy as np
def main():
x = np.array([1,0,0])
y = np.array([0,1,0])
z = np.array([0,0,1])
axis_in =input('Enter the coordinates of rotation axis:\n')
axis  = np.array([np.float(n)for n in axis_in.split()])
#print (axis)
angle_in =input('Enter the rotation angle:\n')
angle  = np.float(angle_in)*np.pi/180
#print (np.degrees(angle))
point_in =input('Enter the coordinates to be rotated:\n')
point = np.array([np.float(n)for n in point_in.split()])
#print (point)
axis1 = np.array([axis[0],axis[1],0])
theta = np.arccos((axis1@x)/(axis1))
print('theta: ', np.degrees(theta))
t1 = np.array([[np.cos(-theta),-np.sin(-theta),0],[np.sin(-theta), np.cos(-theta),0],[0,0,1]])#z print('coordinates transform matrix:\n', t1)
point1 = t1@point
print('transformed coordinates I:\n', point1)
phi = np.arccos((axis@z)/(axis))
print('phi: ', np.degrees(phi))
t2 = np.array([[np.cos(-phi),0, np.sin(-phi)],[0,1,0],[-np.sin(-phi),0, np.cos(-phi)]])#y
print('coordinates transform matrix:\n', t2)
point2 = t2@point1
print('transformed coordinates II:\n', point2)
t3 = np.array([[np.cos(angle),-np.sin(angle),0],[np.sin(angle), np.cos(angle),0],[0,0,1]])#z print('rotation matrix under new coordinates:\n', t3)
rotation = np.transpose(t1) @ np.transpose(t2) @ t3 @ t2 @ t1
print('rotation matrix under original coordinates:\n', rotation)
point_rotated = rotation @ np.transpose(point)
print('the coordinates of rotated point is: \n', point_rotated)
if __name__ =="__main__":
main()
python绘制球⾯渐开线
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
r_theta=np.radians(45)
l=1000
cr,id[0:l*np.sin(r_theta):1000j,0:6.28:50j]
cx=s(ctheta)
cy=cr*np.sin(ctheta)
cz=np.sqrt(cx**2+cy**2)/np.tan(r_theta)
gammar=np.linspace(0,1,100)
theta=np.arctan(np.sin(r_theta)*np.tan(gammar))/np.sin(r_theta)-gammar R=np.arctan(np.tan(r_theta)/np.cos(gammar ))
x=l*np.sin(R)*np.sin(theta)
y=l*np.sin(R)*np.cos(theta)
z=s(R)
fig=plt.figure()
fig.suptitle('sperical Involute')
ax = Axes3D(fig)
ax.plot_surface(cx,cy,cz)
ax.plot(x,y,z)
plt.axis('equal')
plt.show()
多元酸⽔解pH值计算
磷酸⼆氢钾的⽔解系数
化简得到的平衡⽅程：
[H+]^5
+ (k1+c) [H+]^4
+ (k1k2-k) [H+]^3
+ (k1k2k3-(kk1+ck1k2)) [H+]^2
- (kk1k2+2ck1k2k3) [H+]^1
- kk1k2k3 = 0
python求解，有点⽑病，python算不出来，后来改⽤的MATLAB。from sympy import*
x = Symbol('x')
k1 =0.00752
k2 =0.0000000623
k3 =0.00000000000022
k =0.00000000000001
c =0.01
p4=k1+c
linspace函数pythonp3=k1*k2-k
p2=k1*k2*k3-(k*k1+c*k1*k2)
p1=-(k*k1*k2+2*c*k1*k2*k3)
p0=-k*k1*k2*k3
print(p4,p3,p2,p1,p0)
print(solve([x**5+p4*x**4+p3*x**3+p2*x**2+ p1*x+0.1+p0],[x]))

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。