基于Python3的科学运算与常用算法-第3章求解方程

最新推荐文章于 2024-05-04 04:17:15 发布

PhyiscsLogical

最新推荐文章于 2024-05-04 04:17:15 发布

阅读量513

点赞数 3

分类专栏： scipy python numpy

本文链接：https://blog.csdn.net/phyiscslogical/article/details/101515609

版权

python 同时被 3 个专栏收录

12 篇文章 0 订阅

订阅专栏

numpy

8 篇文章 2 订阅

订阅专栏

scipy

5 篇文章 0 订阅

订阅专栏

本系列是使用Ubuntu 64位操作系统上，anaconda环境下，采用jupyter notebook 实现python3科学计算的文章，应用于青岛科技大学信息学院科学计算与数据挖掘等多个课程中python教学。需要安装的模块，numpy ，matplotlib，scipy，sympy，pandas等科学计算常用模块。由于匆忙成稿，个人能力所限，显得粗糙，并且有几个地方是错误的。主要参考了Numerical Python 2nd Edition
————————————————

测试scipy

%time
import numpy as np
from scipy.linalg import solve
#输出系数矩阵
a=np.array([[3,1,-2],[1,-1,4],[2,0,3]])
#值
b=np.array([5,-2,2.5])
#计算
x=solve(a,b)
#打印结果
print(x)

CPU times: user 2 µs, sys: 1 µs, total: 3 µs
Wall time: 3.81 µs
[0.5 4.5 0.5]

测试sympy

from sympy import *
x,y= symbols('x,y')
print(x + y)

x + y

from sympy import *
x= symbols('x')
print(solve(x*2-4,x))

[2]

from sympy import *
x,y= symbols('x,y')
print(solve([2*x-y-3,3*x+y-7],[x,y]))

{x: 2, y: 1}

from sympy import *
x,a= symbols('x,a')
#计算多项式的结果,表达式方式的结果
result_or=solve(x**2+a**2,x)
print(result_or)
#带入数值,这里设置a=3
for temp in result_or:
    print(temp.subs([(a,3)]))

[-I*a, I*a]
-3*I
3*I

ODE

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import sympy

sympy.init_printing()  #为了更好地显示来自 SymPy 的输出，我们需要初始化其屏幕打印系统：

t, k, T0, Ta = sympy.symbols('t, k, T_0, T_a')

T = sympy.Function('T')

ode = T(t).diff(t) + k*(T(t) - Ta)

sympy.Eq(ode)

$\displaystyle k \left(- T_{a} + T{\left(t \right)}\right) + \frac{d}{d t} T{\left(t \right)} = 0$

ode_sol = sympy.dsolve(ode)

ode_sol

$\displaystyle T{\left(t \right)} = C_{1} e^{- k t} + T_{a}$

$\frac{d^2x(t)}{dx^2} + 2\gamma \omega_0 \frac{dx(t)}{dt} + \omega_0^2x(t) = 0$

t, omega0, gamma = sympy.symbols("t, omega_0, gamma", positive=True)
x = sympy.Function('x')
ode = x(t).diff(t, 2) + 2 * gamma * omega0 * x(t).diff(t) + omega0**2*x(t)

sympy.Eq(ode)

$\displaystyle 2 \gamma \omega_{0} \frac{d}{d t} x{\left(t \right)} + \omega_{0}^{2} x{\left(t \right)} + \frac{d^{2}}{d t^{2}} x{\left(t \right)} = 0$

ode_sol = sympy.dsolve(ode)

ode_sol

$\displaystyle x{\left(t \right)} = C_{1} e^{\omega_{0} t \left(- \gamma - \sqrt{\gamma - 1} \sqrt{\gamma + 1}\right)} + C_{2} e^{\omega_{0} t \left(- \gamma + \sqrt{\gamma - 1} \sqrt{\gamma + 1}\right)}$

初值问题

由于这是一个二阶 ODE，因此在通用解中有两个待定积分常数。我们需要为位置$ x(0)$ 和速度 $\frac{dx(t)}{dt}\Bigr|_{t=0}$ 指定初始条件，以便为 ODE 指定一个特解。为此，创建一个包含这些初始条件的字典，并使用 apply_ics 将其应用于 ODE 解：

ics = {x(0): 1, x(t).diff(t).subs(t, 0): 0}

ics

$KaTeX parse error: Undefined control sequence: \substack at position 99: …ght)} \right|_{\̲s̲u̲b̲s̲t̲a̲c̲k̲{ t=0 }} : 0\ri…$

def apply_ics(sol, ics, x, known_params):
    """
    Apply the initial conditions (ics), given as a dictionary on
    the form ics = {y(0): y0, y(x).diff(x).subs(x, 0): yp0, ...},
    to the solution of the ODE with independent variable x.
    The undetermined integration constants C1, C2, ... are extracted
    from the free symbols of the ODE solution, excluding symbols in
    the known_params list.
    """

    free_params = sol.free_symbols - set(known_params)
    eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics)
            for n in range(len(ics))]
    sol_params = sympy.solve(eqs, free_params)

    return sol.subs(sol_params)

 x_t_sol =  apply_ics(ode_sol, ics, t, [omega0, gamma])

x_t_sol

$\displaystyle x{\left(t \right)} = \left(- \frac{\gamma}{2 \sqrt{\gamma^{2} - 1}} + \frac{1}{2}\right) e^{\omega_{0} t \left(- \gamma - \sqrt{\gamma - 1} \sqrt{\gamma + 1}\right)} + \left(\frac{\gamma}{2 \sqrt{\gamma^{2} - 1}} + \frac{1}{2}\right) e^{\omega_{0} t \left(- \gamma + \sqrt{\gamma - 1} \sqrt{\gamma + 1}\right)}$

对于 $\gamma = 1$ ，通常计算会产生零除法错误，所以处理的方法是在这个值范围内要，采用极限。

x_t_critical = sympy.limit(x_t_sol.rhs, gamma, 1)

x_t_critical

$\displaystyle \left(\omega_{0} t + 1\right) e^{- \omega_{0} t}$

最后，我们绘制 ω0=2π 和阻尼比 γ 不同值的图形：

fig, ax = plt.subplots(figsize=(8, 4))
tt = np.linspace(0, 3, 250)
w0 = 2 * sympy.pi
for g in [0.1, 0.5, 1, 2.0, 5.0]:
    if g == 1:
        x_t = sympy.lambdify(t, x_t_critical.subs({omega0: w0}), 'numpy')
    else:
        x_t = sympy.lambdify(t, x_t_sol.rhs.subs({omega0: w0, gamma: g}), 'numpy')
    ax.plot(tt, x_t(tt).real, label=r"$\gamma = %.1f$" % g)
ax.set_xlabel(r"$t$", fontsize=18)
ax.set_ylabel(r"$x(t)$", fontsize=18)
ax.legend()

<matplotlib.legend.Legend at 0x7fcfe0caea50>

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-n4ccT5Tf-1569550074917)(output_34_1.png)]

参考文献

https://blog.csdn.net/and_w/article/details/80160033
https://blog.csdn.net/And_w/article/details/80160033
https://www.cnblogs.com/xzcfightingup/p/7598293.html
https://blog.csdn.net/meiqi0538/article/details/82990432
https://www.jianshu.com/p/339c91ae9f41
https://www.cnblogs.com/wj-1314/p/10244807.html
https://www.jianshu.com/p/b2abedd3784e
https://www.jianshu.com/p/339c91ae9f41