Python科学计算入门

alt text

1. 为什么选择Python？

Python已成为科学计算和数据分析的首选语言，原因包括：

  • 语法简洁易学
  • 丰富的科学计算库
  • 强大的可视化能力
  • 活跃的社区支持

2. 核心库介绍

2.1 NumPy - 数值计算基础

NumPy提供高性能的多维数组对象和相关工具。

import numpy as np

# 创建数组
arr = np.array([1, 2, 3, 4, 5])

# 数组运算
arr_squared = arr ** 2

# 矩阵运算
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = np.dot(A, B)  # 矩阵乘法

2.2 SciPy - 科学计算工具

SciPy基于NumPy，提供更高级的科学计算功能：

from scipy import optimize, integrate

# 求解非线性方程
def func(x):
    return x**2 - 4

root = optimize.fsolve(func, 1)

# 数值积分
result, error = integrate.quad(lambda x: x**2, 0, 1)

2.3 Matplotlib - 数据可视化

import matplotlib.pyplot as plt

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

plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('sin(x)')
plt.title('正弦函数')
plt.grid(True)
plt.show()

3. 工程应用实例

3.1 应力分析

import numpy as np

# 材料参数
E = 200e9  # 弹性模量 (Pa)
nu = 0.3   # 泊松比

# 应变
epsilon_x = 0.001
epsilon_y = 0.0005

# 应力计算（平面应力）
sigma_x = E / (1 - nu**2) * (epsilon_x + nu * epsilon_y)
sigma_y = E / (1 - nu**2) * (epsilon_y + nu * epsilon_x)

print(f"σx = {sigma_x/1e6:.2f} MPa")
print(f"σy = {sigma_y/1e6:.2f} MPa")

3.2 梁的挠度计算

import numpy as np
import matplotlib.pyplot as plt

# 简支梁中点集中载荷
L = 5.0    # 梁长度 (m)
P = 10000  # 载荷 (N)
E = 200e9  # 弹性模量 (Pa)
I = 1e-4   # 惯性矩 (m^4)

# 挠度方程
def deflection(x):
    if x <= L/2:
        return -P*x/(48*E*I) * (3*L**2 - 4*x**2)
    else:
        return -P*(L-x)/(48*E*I) * (3*L**2 - 4*(L-x)**2)

# 计算并绘图
x = np.linspace(0, L, 100)
y = [deflection(xi) for xi in x]

plt.plot(x, y)
plt.xlabel('位置 (m)')
plt.ylabel('挠度 (m)')
plt.title('简支梁挠度曲线')
plt.grid(True)
plt.show()

4. 数值求解微分方程

4.1 一阶ODE

使用scipy.integrate.odeint求解：

from scipy.integrate import odeint

# 定义微分方程 dy/dt = -k*y
def model(y, t, k):
    dydt = -k * y
    return dydt

# 初始条件
y0 = 5
k = 0.3
t = np.linspace(0, 20, 100)

# 求解
y = odeint(model, y0, t, args=(k,))

# 绘图
plt.plot(t, y)
plt.xlabel('时间')
plt.ylabel('y(t)')
plt.title('一阶ODE数值解')
plt.show()

5. 最优化问题

5.1 无约束优化

from scipy.optimize import minimize

# 目标函数
def objective(x):
    return x[0]**2 + x[1]**2 + x[0]*x[1]

# 初始猜测
x0 = [1, 1]

# 优化
result = minimize(objective, x0, method='BFGS')
print(f"最优解: x = {result.x}")
print(f"最小值: f(x) = {result.fun}")

5.2 约束优化

from scipy.optimize import minimize

def objective(x):
    return -x[0] * x[1]  # 最大化 x*y

def constraint1(x):
    return x[0] + 2*x[1] - 10  # x + 2y <= 10

con1 = {'type': 'ineq', 'fun': constraint1}
bounds = [(0, None), (0, None)]

result = minimize(objective, [1, 1],
                  method='SLSQP',
                  bounds=bounds,
                  constraints=con1)

6. 符号计算 - SymPy

对于需要符号推导的问题：

from sympy import symbols, diff, integrate, simplify

x, y = symbols('x y')

# 求导
f = x**3 + 2*x**2 - x
df = diff(f, x)
print(f"f'(x) = {df}")

# 积分
integral = integrate(f, x)
print(f"∫f(x)dx = {integral}")

# 化简
expr = (x**2 - 1)/(x - 1)
simplified = simplify(expr)
print(f"化简: {simplified}")

7. 推荐学习路径

  1. 基础：Python语法 → NumPy数组操作
  2. 进阶：SciPy科学计算 → Matplotlib可视化
  3. 应用：工程问题建模 → 数值求解
  4. 扩展：Pandas数据分析 → SymPy符号计算

8. 常用资源

  • 官方文档：numpy.org, scipy.org
  • 教程网站：Real Python, Python for Engineers
  • 书籍：《Python科学计算》、《Numerical Python》

9. 总结

Python为工程师提供了强大而灵活的科学计算工具，掌握这些库可以大大提高工作效率。建议从简单问题入手，逐步深入复杂应用。

下期将介绍如何用Python进行有限元后处理。