跳转至

SciPy

约 3554 个字 386 行代码 预计阅读时间 17 分钟

1

完成示例

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

m1 = 4866
m2 = 2433
k1 = 31557.298
k2 = 80000
k3 = 10000
k4 = 656.3516
F = 6250
w = 1.4005
m_c = 1335.535


# 定义二元微分方程函数
def f(y, x, m1, m2, k1, k2, k3, k4, F, w, m_c):
    # return [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * (y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c),
    #         (k2 * (y[0] - y[1]) + k3 * (y[2] - y[3]) ** 1.5) / m2]
    sign = np.sign(y[2] - y[3])
    return [y[2], y[3],
            (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * sign * abs(y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (
                        m1 + m_c),
            (k2 * (y[0] - y[1]) + k3 * sign * abs(y[2] - y[3]) ** 1.5) / m2]

# 初始条件
y0 = [0, 0, 0, 0]

# 自变量范围
x = np.linspace(0, 5, 1000)

# 求解微分方程组
sol = odeint(f, y0, x, args=(m1, m2, k1, k2, k3, k4, F, w, m_c))

# 分别提取y1和y2的值
y1 = sol[:, 0]
y2 = sol[:, 1]
y3 = sol[:, 2]
y4 = sol[:, 3]

# 绘制函数图像
plt.plot(x, y1, label='x1')
plt.plot(x, y2, label='x2')
plt.plot(x, y3, label='x3')
plt.plot(x, y4, label='x4')

# 添加图例、坐标轴标签和标题
plt.legend()
plt.xlabel('t')
plt.ylabel('x')
plt.title('Solution of a System of Differential Equations')

# 显示图像
plt.show()

2

np.linspace(0, 5, 1000) 

@array_function_dispatch(_linspace_dispatcher)
def linspace(start, stop, num=50, ...):
    ...

第一、二个参数是起始和结束，第三个参数是分割点的数量（越多点越密，曲线越平滑），默认50

chat log

user:

你知道scipy吗？

chatgpt:

是的，我知道Scipy。Scipy是一个用于数学、科学和工程计算的Python库。它建立在NumPy库的基础上，提供了许多用于信号处理、优化、统计分析、线性代数和数值积分等方面的函数和工具。Scipy还包括了一些特殊函数、傅里叶变换、信号滤波器、最优化器、线性代数工具和随机数生成器等功能。Scipy是Python科学计算库的一个强大而广泛使用的组件，其灵活性和易用性使其成为科学家、工程师和研究人员的首选之一。

user:

怎么使用scipy你知道吗？

chatgpt:

是的，我可以给您提供使用Scipy的基本步骤：

  1. 安装Scipy库：在您的Python环境中使用pip命令来安装Scipy库。例如在命令提示符中输入：pip install scipy。

  2. 导入Scipy模块：在Python代码中导入Scipy模块。例如：import scipy。

  3. 使用Scipy函数：使用Scipy提供的函数来进行数学、科学和工程计算。例如，要进行线性代数计算，可以使用Scipy中的linalg模块，如下所示：

import scipy.linalg as la

# 创建一个2x2的矩阵
A = np.array([[1, 2], [3, 4]])

# 计算矩阵的行列式
det_A = la.det(A)

# 计算矩阵的逆
inv_A = la.inv(A)

除了线性代数，Scipy还提供了许多其他模块和函数，例如：

  • scipy.optimize：提供了优化算法，用于寻找函数的最小值或最大值。
  • scipy.integrate：提供了数值积分函数。
  • scipy.stats：提供了统计分析函数，用于生成和操作概率分布。
  • scipy.signal：提供了信号处理函数，例如滤波器和傅里叶变换等。

Scipy库的使用非常广泛，如果您需要进行特定的数学、科学或工程计算，可以先看看Scipy库是否提供了相关的函数和工具。

user:

那你知道如何使用scipy解二元微分方程组吗？

chatgpt:

是的，使用Scipy可以解二元微分方程组。下面是一个简单的例子，展示了如何使用Scipy解二元微分方程组：

假设我们要解决以下二元微分方程组：

y1' = -2y1 + y2 y2' = y1 - 2y2

其中y1和y2是关于自变量x的函数。为了解决这个方程组，我们需要将它转换为矩阵形式：

[d/dx(y1)] [-2 1] [y1] [d/dx(y2)] = [ 1 -2] * [y2]

然后，我们可以使用Scipy的odeint函数来解决该方程组。odeint函数可以用来解决一般形式的常微分方程组。以下是一个使用odeint函数来解决上述二元微分方程组的示例代码：

import numpy as np
from scipy.integrate import odeint

# 定义二元微分方程函数
def f(y, x):
    return [-2*y[0] + y[1], y[0] - 2*y[1]]

# 初始条件
y0 = [1, 0]

# 自变量范围
x = np.linspace(0, 10, 100)

# 求解微分方程组
sol = odeint(f, y0, x)

# 分别提取y1和y2的值
y1 = sol[:, 0]
y2 = sol[:, 1]

在这个示例中，我们首先定义了一个函数f来表示二元微分方程组。然后，我们指定了初始条件y0和自变量范围x，并使用odeint函数来求解微分方程组。最后，我们提取了y1和y2的值，并将它们存储在y1和y2数组中。

请注意，odeint函数的第一个参数是表示微分方程组的函数名，第二个参数是初始条件，第三个参数是自变量范围。在这个例子中，我们使用numpy.linspace函数来生成自变量范围x。

user:

如何打印y1 y2成函数图像？

chatgpt:

要将y1和y2打印成函数图像，我们需要使用Matplotlib库。Matplotlib是一个Python库，可用于绘制各种类型的图形，包括线图、散点图、条形图、饼图等等。

以下是一个示例代码，展示了如何使用Matplotlib库绘制y1和y2的函数图像：

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

# 定义二元微分方程函数
def f(y, x):
    return [-2*y[0] + y[1], y[0] - 2*y[1]]

# 初始条件
y0 = [1, 0]

# 自变量范围
x = np.linspace(0, 10, 100)

# 求解微分方程组
sol = odeint(f, y0, x)

# 分别提取y1和y2的值
y1 = sol[:, 0]
y2 = sol[:, 1]

# 绘制函数图像
plt.plot(x, y1, label='y1')
plt.plot(x, y2, label='y2')

# 添加图例、坐标轴标签和标题
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of a System of Differential Equations')

# 显示图像
plt.show()

在这个示例中，我们首先求解了微分方程组，并提取了y1和y2的值。然后，我们使用Matplotlib的plot函数绘制了y1和y2的函数图像，并使用legend函数添加了图例。最后，我们使用xlabel、ylabel和title函数添加了坐标轴标签和标题，并使用show函数显示了图像。

这个代码会生成一个包含y1和y2函数图像的窗口。您可以通过调整x范围、初始条件和微分方程函数来生成不同的图像。

user:

系数浮点数时，报错了， TypeError: unsupported operand type(s) for *: 'function' and 'float' 应该如何解决？

chatgpt:

这个错误通常是由于参数类型不匹配引起的。具体地说，在您的代码中，您可能在使用odeint函数时将一个函数对象与一个浮点数相乘，从而导致了这个错误。

要解决这个问题，您可以检查您的代码，确保您没有在使用odeint函数时将函数对象与浮点数相乘。另外，您可以尝试将所有的系数都显式地转换为浮点数，以确保odeint函数接收到的参数类型都是一致的。

以下是一个修复这个问题的示例代码：

import numpy as np
from scipy.integrate import odeint

# 定义二元微分方程函数
def f(y, x, a, b, c, d):
    return [-a*y[0] + b*y[1], c*y[0] - d*y[1]]

# 初始条件
y0 = [1, 0]

# 自变量范围
x = np.linspace(0, 10, 100)

# 定义参数
a = 1.0
b = 2.0
c = 3.0
d = 4.0

# 求解微分方程组
sol = odeint(f, y0, x, args=(a, b, c, d))

# 分别提取y1和y2的值
y1 = sol[:, 0]
y2 = sol[:, 1]

在这个示例中，我们首先将所有的系数都显式地转换为浮点数。然后，我们将这些系数作为参数传递给二元微分方程函数f，并在调用odeint函数时将这些参数作为附加参数传递。最后，我们求解微分方程组，并提取了y1和y2的值。

请注意，我们将参数列表(a, b, c, d)作为odeint函数的args参数传递。这是因为，odeint函数需要接收一个函数和一个初始条件，并可选地接收一些附加的参数。这些附加参数必须作为一个元组传递给args参数，以便在求解微分方程组时使用。

user:

那么scipy如何求解三元二阶微分方程？

chatgpt:

Scipy库中没有直接求解三元二阶微分方程的函数，但可以使用odeint函数来求解。为了使用odeint函数求解三元二阶微分方程，需要将其转化为一个六元一阶微分方程组。

例如，假设我们要求解以下三元二阶微分方程组：

x''(t) + 2y'(t) + 3z(t) = 0 y''(t) + 3x(t) + 4z'(t) = 0 z''(t) + x(t) + 2*y(t) = 0

我们可以将它转化为以下六元一阶微分方程组：

x'(t) = v_x(t) v_x'(t) = -2v_y(t) - 3z(t) y'(t) = v_y(t) v_y'(t) = -3x(t) - 4v_z(t) z'(t) = v_z(t) v_z'(t) = -x(t) - 2*y(t)

其中，x(t)、y(t)和z(t)是待求解的函数，v_x(t)、v_y(t)和v_z(t)是它们的导数。

我们可以将这个六元一阶微分方程组作为odeint函数的输入，并使用初始条件求解它。以下是一个使用odeint函数求解上述三元二阶微分方程的示例代码：

import numpy as np
from scipy.integrate import odeint

# 定义六元一阶微分方程函数
def f(y, t):
    x, vx, y, vy, z, vz = y
    return [vx, -2*vy - 3*z, vy, -3*x - 4*vz, vz, -x - 2*y]

# 初始条件
y0 = [1, 0, 1, 0, 1, 0]

# 自变量范围
t = np.linspace(0, 10, 100)

# 求解微分方程组
sol = odeint(f, y0, t)

# 分别提取x, y, z 的值
x = sol[:, 0]
y = sol[:, 2]
z = sol[:, 4]

在这个示例中，我们定义了一个函数f来表示六元一阶微分方程组，并指定了初始条件y0和自变量范围t。然后，我们使用odeint函数求解微分方程组，并提取了x、y和z的值。

请注意，在这个例子中，我们将初始条件y0作为一个六元列表传递给odeint函数。在函数f中，我们使用y列表来表示六个函数和它们的导数。

user:

如何打印成图像？

chatgpt:

要将求解的三元二阶微分方程的结果打印成图像，可以使用Matplotlib库。Matplotlib是一个Python库，可用于绘制各种类型的图形，包括线图、散点图、条形图、饼图等等。

以下是一个示例代码，展示了如何使用Matplotlib库绘制x、y和z的函数图像：

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

# 定义六元一阶微分方程函数
def f(y, t):
    x, vx, y, vy, z, vz = y
    return [vx, -2*vy - 3*z, vy, -3*x - 4*vz, vz, -x - 2*y]

# 初始条件
y0 = [1, 0, 1, 0, 1, 0]

# 自变量范围
t = np.linspace(0, 10, 100)

# 求解微分方程组
sol = odeint(f, y0, t)

# 分别提取x, y, z的值
x = sol[:, 0]
y = sol[:, 2]
z = sol[:, 4]

# 绘制函数图像
plt.plot(t, x, label='x')
plt.plot(t, y, label='y')
plt.plot(t, z, label='z')

# 添加图例、坐标轴标签和标题
plt.legend()
plt.xlabel('t')
plt.ylabel('f(t)')
plt.title('Solution of a System of Differential Equations')

# 显示图像
plt.show()

在这个示例中，我们首先求解了微分方程组，并提取了x、y和z的值。然后，我们使用Matplotlib的plot函数绘制了x、y和z的函数图像，并使用legend函数添加了图例。最后，我们使用xlabel、ylabel和title函数添加了坐标轴标签和标题，并使用show函数显示了图像。

这个代码会生成一个包含x、y和z函数图像的窗口。您可以通过调整自变量范围、初始条件和微分方程函数来生成不同的图像。

user:

能不能打印成3d的图像？

chatgpt:

是的，可以使用Matplotlib库中的mplot3d模块将三元二阶微分方程的求解结果打印成3D图像。

以下是一个示例代码，展示了如何使用Matplotlib的mplot3d模块绘制x、y和z的3D图像：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 定义六元一阶微分方程函数
def f(y, t):
    x, vx, y, vy, z, vz = y
    return [vx, -2*vy - 3*z, vy, -3*x - 4*vz, vz, -x - 2*y]

# 初始条件
y0 = [1, 0, 1, 0, 1, 0]

# 自变量范围
t = np.linspace(0, 10, 100)

# 求解微分方程组
sol = odeint(f, y0, t)

# 分别提取x, y, z的值
x = sol[:, 0]
y = sol[:, 2]
z = sol[:, 4]

# 绘制3D图像
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z)

# 添加坐标轴标签和标题
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.title('Solution of a System of Differential Equations')

# 显示图像
plt.show()

在这个示例中，我们首先求解了微分方程组，并提取了x、y和z的值。然后，我们使用Matplotlib的mplot3d模块中的Axes3D类创建一个3D坐标系，并使用plot函数绘制了x、y和z的3D图像。最后，我们使用set_xlabel、set_ylabel和set_zlabel函数添加了坐标轴标签，并使用title函数添加了标题。

这个代码会生成一个包含x、y和z函数的3D图像的窗口。您可以通过调整自变量范围、初始条件和微分方程函数来生成不同的图像。

user:

那scipy能不能解三元非线性微分方程？

chatgpt:

Scipy库中的odeint函数可以用于求解三元非线性微分方程。

对于三元非线性微分方程，可以将其转化为一个六元一阶微分方程组，并将其作为odeint函数的输入。具体来说，如果我们有以下三元非线性微分方程：

x''(t) = -x(t) * y(t) * z(t) y''(t) = -y(t) * z(t) * x(t) z''(t) = -z(t) * x(t) * y(t)

我们可以将其转化为以下六元一阶微分方程组：

x'(t) = v_x(t) v_x'(t) = -x(t) * y(t) * z(t) y'(t) = v_y(t) v_y'(t) = -y(t) * z(t) * x(t) z'(t) = v_z(t) v_z'(t) = -z(t) * x(t) * y(t)

其中，x(t)、y(t)和z(t)是待求解的函数，v_x(t)、v_y(t)和v_z(t)是它们的导数。

以下是一个使用odeint函数求解上述三元非线性微分方程的示例代码：

import numpy as np
from scipy.integrate import odeint

# 定义六元一阶微分方程函数
def f(y, t):
    x, vx, y, vy, z, vz = y
    return [vx, -x*y*z, vy, -y*z*x, vz, -z*x*y]

# 初始条件
y0 = [1, 0, 1, 0, 1, 0]

# 自变量范围
t = np.linspace(0, 10, 100)

# 求解微分方程组
sol = odeint(f, y0, t)

# 分别提取x, y, z 的值
x = sol[:, 0]
y = sol[:, 2]
z = sol[:, 4]

在这个示例中，我们定义了一个函数f来表示六元一阶微分方程组，并指定了初始条件y0和自变量范围t。然后，我们使用odeint函数求解微分方程组，并提取了x、y和z的值。

请注意，在这个例子中，我们将初始条件y0作为一个六元列表传递给odeint函数。在函数f中，我们使用y列表来表示六个函数和它们的导数。

需要注意的是，对于非线性微分方程，求解可能会比较困难，需要选择合适的方法和参数来求解。如果遇到求解困难的情况，可以尝试使用其他数值求解方法，或者手动调整odeint函数的参数来提高求解效果。

user: 

RuntimeWarning: invalid value encountered in scalar power
(k2 * (y[0] - y[1]) + k3 * (y[2] - y[3]) ** 1.5) / m2]
这个报错应该如何解决？ 完整代码如下： 
import numpy as np
from scipy import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from math import *

m1 = 4.866
m2 = 2433
k1 = 31557.298
k2 = 80000
k3 = 10000
k4 = 656.3516
F = 6250
w = 1.4005
m_c = 1335.535


# 定义二元微分方程函数
def f(y, x, m1, m2, k1, k2, k3, k4, F, w, m_c):
    return [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * (y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c),
            (k2 * (y[0] - y[1]) + k3 * (y[2] - y[3]) ** 1.5) / m2]


# 初始条件
y0 = [0, 0, 0, 0]

# 自变量范围
x = np.linspace(0, 5, 1000)

# 求解微分方程组
sol = odeint(f, y0, x, args=(m1, m2, k1, k2, k3, k4, F, w, m_c))

# 分别提取y1和y2的值
y1 = sol[:, 0]
y2 = sol[:, 1]
y3 = sol[:, 2]
y4 = sol[:, 3]

# 绘制函数图像
plt.plot(x, y1, label='y1')
plt.plot(x, y2, label='y2')
plt.plot(x, y3, label='y3')
plt.plot(x, y4, label='y4')

# 添加图例、坐标轴标签和标题
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of a System of Differential Equations')

# 显示图像
plt.show()

chatgpt:

这个报错是因为在求解微分方程组时，出现了非法值。具体来说，在函数f中的第四个方程中，当分母m2为0时，会发生除零错误，进而导致错误的计算结果。

为了解决这个问题，可以在函数f中对分母m2进行判断，如果m2为0，则将y[3]的导数设为0。这样可以避免除零错误，并保证求解结果的正确性。修改后的代码如下：

# 定义二元微分方程函数
def f(y, x, m1, m2, k1, k2, k3, k4, F, w, m_c):
    if m2 == 0:
        dydt = [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * (y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c), 0]
    else:
        dydt = [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * (y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c),
                (k2 * (y[0] - y[1]) + k3 * (y[2] - y[3]) ** 1.5) / m2]
    return dydt

在修改后的代码中，我们首先对m2进行判断，如果m2为0，则将y[3]的导数设为0，否则按原来的方式计算y[3]的导数。这样可以避免除零错误，并保证求解结果的正确性。

需要注意的是，在实际应用中，出现非法值的原因可能不仅仅是除零错误，还可能包括其他情况。因此，在编写程序时，需要对可能出现的非法值进行充分的判断和处理，以保证程序的正确性和稳定性。

user: 

RuntimeWarning: invalid value encountered in scalar power
(k2 * (y[0] - y[1]) + k3 * (y[2] - y[3]) ** 1.5) / m2]
How should this error be resolved? The complete code is as follows: 
import numpy as np
from scipy import *
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from math import *

m1 = 4.866
m2 = 2433
k1 = 31557.298
k2 = 80000
k3 = 10000
k4 = 656.3516
F = 6250
w = 1.4005
m_c = 1335.535


# 定义二元微分方程函数
def f(y, x, m1, m2, k1, k2, k3, k4, F, w, m_c):
    return [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * (y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c),
            (k2 * (y[0] - y[1]) + k3 * (y[2] - y[3]) ** 1.5) / m2]


# 初始条件
y0 = [0, 0, 0, 0]

# 自变量范围
x = np.linspace(0, 5, 1000)

# 求解微分方程组
sol = odeint(f, y0, x, args=(m1, m2, k1, k2, k3, k4, F, w, m_c))

# 分别提取y1和y2的值
y1 = sol[:, 0]
y2 = sol[:, 1]
y3 = sol[:, 2]
y4 = sol[:, 3]

# 绘制函数图像
plt.plot(x, y1, label='y1')
plt.plot(x, y2, label='y2')
plt.plot(x, y3, label='y3')
plt.plot(x, y4, label='y4')

# 添加图例、坐标轴标签和标题
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Solution of a System of Differential Equations')

# 显示图像
plt.show()

chatgpt:

This error occurs because an invalid value was encountered during the calculation of a scalar power operation. Specifically, in the fourth equation of the function f, when (y[2] - y[3]) is negative, the power of 1.5 results in an imaginary number, which is not a valid value for a real scalar. This leads to the RuntimeWarning error.

To fix this error, you can use the numpy.sign function to get the sign of the quantity (y[2] - y[3]), and then multiply it with the absolute value of this quantity raised to the power of 1.5. The resulting expression will always be real and prevent this error from occurring. Here's the modified code:

# Define the system of differential equations
def f(y, x, m1, m2, k1, k2, k3, k4, F, w, m_c):
    sign = np.sign(y[2] - y[3])
    return [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * sign * abs(y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c),
            (k2 * (y[0] - y[1]) + k3 * sign * abs(y[2] - y[3]) ** 1.5) / m2]

In this modified code, we first compute the sign of (y[2] - y[3]) using numpy.sign. We then use this sign to multiply the absolute value of (y[2] - y[3]) raised to the power of 1.5. This ensures that the expression is always real, regardless of the sign of (y[2] - y[3]), and resolves the RuntimeWarning error.

Note that this error can occur in other situations where a scalar power operation results in an invalid value, so it's always a good idea to check for such cases and handle them appropriately to avoid errors.

user:

你能用中文再输出一遍吗？

chatgpt:

这个错误是因为在计算标量乘方时出现了非法值。具体来说，在函数f的第四个方程中，当(y[2] - y[3])为负数时，将1.5次幂应用于该负数会导致虚数，这不是标量的有效值，因此会出现RuntimeWarning错误。

为了解决这个问题，可以使用numpy.sign函数获取(y[2] - y[3])的符号，并将其与绝对值的1.5次幂相乘，可以确保计算结果总是实数，从而避免出现错误。以下是修改后的代码：

# 定义二元微分方程函数
def f(y, x, m1, m2, k1, k2, k3, k4, F, w, m_c):
    sign = np.sign(y[2] - y[3])
    return [y[2], y[3], (-k1 * y[0] - k2 * (y[0] - y[1]) - k3 * sign * abs(y[2] - y[3]) ** 1.5 - k4 * y[2] + F * np.cos(w * x)) / (m1 + m_c),
            (k2 * (y[0] - y[1]) + k3 * sign * abs(y[2] - y[3]) ** 1.5) / m2]

在修改后的代码中，我们首先使用numpy.sign获取(y[2] - y[3])的符号，然后将其与绝对值的1.5次幂相乘。这将确保计算结果总是实数，并避免出现RuntimeWarning错误。

需要注意的是，这种错误可能会在其他情况下出现，因此，在编写程序时，需要检查这些情况并进行适当的处理，以避免出现错误。

最后更新: 2023-10-24
创建日期: 2023-07-17