常数变异法求二阶微分方程

### 使用常数变异法解二阶微分方程

对于二阶线性非齐次微分方程 \( y'' + p(x)y' + q(x)y = f(x) \)，可以通过常数变异法来求解。该方法基于已知对应的齐次方程 \( y'' + p(x)y' + q(x)y = 0 \) 的通解。

假设齐次方程有两个线性无关的特解 \( y_1(x) \) 和 \( y_2(x) \)，那么其通解形式为：

\[ Y_h(x) = c_1y_1(x) + c_2y_2(x) \]

为了找到非齐次方程的一个特解，设未知函数的形式如下：

\[ Y_p(x) = u_1(x)y_1(x) + u_2(x)y_2(x) \]

其中 \( u_1(x) \) 和 \( u_2(x) \) 是待定的新变量。接下来需要建立关于这两个新的未知函数的一组辅助条件以便简化问题。通常采用拉格朗日公式作为附加约束:

\[ u'_1(x)y_1(x) + u'_2(x)y_2(x)=0 \]

这样做的目的是让导数项只保留一阶导数部分，从而使得原方程降阶处理成为可能。将上面设定带入原始非齐次方程并利用给定关系消去高阶导数之后可得另一个表达式涉及\(u'\):

\[ u'_1(x)y'_1(x)+u'_2(x)y'_2(x)=f(x)/W(y_1,y_2)(x)\]

这里引入了 Wronskian 行列式的概念 \( W(y_1, y_2)(x) = y_1(x)y'_2(x)-y'_1(x)y_2(x) \)[^1]。

联立以上两式构成一个简单的线性系统，可以直接积分求出 \( u_1(x), u_2(x) \)，进而获得完整的非齐次方程的通解。

from sympy import symbols, Function, Eq, dsolve, diff, integrate

# 定义符号和函数
x = symbols('x')
y = Function('y')(x)

p, q, r = map(Function, 'pqr') # 方程系数
c1, c2 = symbols('c1 c2')

# 已知齐次解
Yh = c1 * exp(-integrate(p(x)/2,x)) + c2*x*exp(-integrate(p(x)/2,x))

# 特殊情况下选取的基础解系 (具体取决于实际问题中的p,q,r)
y1 = exp(-integrate(p(x)/2,x))
y2 = x*y1

# 计算wroskian行列式
Wr = y1.diff(x)*y2 - y1*y2.diff(x)

# 构造变易后的方程组
eq_u1_prime = Eq(diff(c1*(x)), (-r(x)*y2/Wr).simplify())
eq_u2_prime = Eq(diff(c2*(x)), (r(x)*y1/Wr).simplify())

# 积分得到u1,u2的具体表示
sol_c1 = integrate((-r(x)*y2/Wr).simplify(), x)
sol_c2 = integrate((r(x)*y1/Wr).simplify(), x)

print(f'u1\'(x): {eq_u1_prime}')
print(f'u2\'(x): {eq_u2_prime}')

# 非齐次特解
yparticular = sol_c1.subs(c1(sol_c1.free_symbols.pop()),symbols('C1'))*y1 \
            + sol_c2.subs(c2(sol_c2.free_symbols.pop()),symbols('C2'))*y2

full_solution = Yh + yparticular

此代码片段展示了如何使用 Python SymPy 库实现上述理论过程的一部分自动化计算。注意这只是一个示意性的例子，在真实的应用场景下还需要根据具体的边界条件调整初始猜测值以及其他细节。