用Python运算常系数齐次线性微分方程的代码

可以使用SymPy库来求解常系数齐次线性微分方程。

首先，我们需要导入SymPy库：

import sympy as sp

然后，定义未知函数和自变量：

x = sp.Symbol('x')
y = sp.Function('y')(x)

接下来，定义常系数齐次线性微分方程：

eq = sp.Eq(y.diff(x, 2) + 2*y.diff(x) + y, 0)

其中，`y.diff(x, 2)` 表示对 `y` 求二阶导数，`y.diff(x)` 表示对 `y` 求一阶导数。

最后，使用 `dsolve` 函数求解微分方程：

sol = sp.dsolve(eq, y)

完整代码如下：

import sympy as sp

x = sp.Symbol('x')
y = sp.Function('y')(x)

eq = sp.Eq(y.diff(x, 2) + 2*y.diff(x) + y, 0)

sol = sp.dsolve(eq, y)

print(sol)

输出结果为：

Eq(y(x), (C1 + C2*x)*exp(-x))

其中，`C1` 和 `C2` 是任意常数，代表微分方程的通解。

相关问题

不使用第三方库用Python运算常系数齐次线性微分方程的代码

要解决常系数齐次线性微分方程，可以使用以下代码：

def solve_homogeneous_linear_ode(a, b, c):
    """
    解常系数齐次线性微分方程 a*y'' + b*y' + c*y = 0
    返回形如 C1*e^(r1*x) + C2*e^(r2*x) 的通解，其中 r1, r2 为方程的两个解
    """
    # 计算判别式
    D = b**2 - 4*a*c
    if D > 0:
        # 有两个不相等的实根
        r1 = (-b + math.sqrt(D)) / (2*a)
        r2 = (-b - math.sqrt(D)) / (2*a)
        return lambda x: C1*math.exp(r1*x) + C2*math.exp(r2*x)
    elif D == 0:
        # 有一个重根
        r = -b / (2*a)
        return lambda x: (C1 + C2*x)*math.exp(r*x)
    else:
        # 有两个共轭的复根
        alpha = -b / (2*a)
        beta = math.sqrt(-D) / (2*a)
        return lambda x: math.exp(alpha*x)*(C1*math.cos(beta*x) + C2*math.sin(beta*x))

其中，`a`、`b`、`c` 分别为微分方程中的系数，函数返回的是一个 lambda 表达式，它可以计算出微分方程的通解。

需要注意的是，上面的代码只能解决常系数齐次线性微分方程，如果方程中有非常数项或不是齐次的，需要另写代码处理。

python 微分方程ipynb

### Python 解决微分方程的 Jupyter Notebook 示例

#### 使用 `SciPy` 库求解常微分方程 (ODE)

为了展示如何利用 Python 和 Jupyter Notebook 来解决微分方程，下面给出一个简单的例子——洛伦兹吸引子模型。该模型描述的是流体动力学中的一种混沌现象。

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def lorenz(t, xyz, sigma=10., beta=8./3, rho=28.0):
    """The Lorenz equations."""
    x, y, z = xyz
    return [sigma * (y - x),
            x * (rho - z) - y,
            x * y - beta * z]

# Initial state at time t=0 is (-8, 8, 27)
initial_state = [-8, 8, 27]
time_span = [0, 40]  # Time interval to integrate over.
solution = solve_ivp(lorenz, time_span, initial_state, dense_output=True)

times = np.linspace(time_span[0], time_span[-1], 10000)
xyz_sol = solution.sol(times)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot(xyz_sol[0], xyz_sol[1], xyz_sol[2])
plt.show()

这段代码展示了怎样定义并绘制三维空间内的轨迹图来表示系统的演化过程[^1]。

#### 利用 `SymPy` 进行符号运算解析解法

除了数值方法外，还可以借助 SymPy 工具包来进行代数操作获得某些类型的 ODE 的精确表达式：

from sympy import *
init_printing(use_unicode=True)

t = symbols('t')  # Define symbol 't'
x = Function('x')(t)  # Create function object named "x"

ode_example = Eq(diff(x,t), sin(t))  # Example of simple first order linear differential equation dx/dt=sin(t).
dsolve(ode_example,x).doit()         # Solve the given ordinary differential equation and display result.

上述脚本会返回给定一阶线性齐次微分方程 \( \frac{dx}{dt}=\sin{(t)}\) 的通解形式[^2]。