sympy dsolve returns 错误答案
sympy dsolve returns incorrect answer
我正在使用 sympy.dsolve 求解衰减链的简单 ODE。
对于不同的衰减率(例如 lambda_1 > lambda_2),我得到的答案是错误的。代入 C1=0 后,我得到一个简单的指数
-N_0*lambda_1*exp(-lambda_1*t)/(lambda_1 - lambda_2)
而不是正确答案有:
(exp(-lambda_1*t)-exp(-lambda_2*t)).
我做错了什么?
这是我的代码
sp.var('lambda_1,lambda_2 t')
sp.var('N_0')
daughter = sp.Function('N2',Positive=True)(t)
stage1 = N_0*sp.exp(-lambda_1*t)
eq = sp.Eq(daughter.diff(t),stage1*lambda_1 - daughter*lambda_2)
sp.dsolve(eq,daughter)
您的微分方程是(使用较短的变量标识符)
y' = A*N0*exp(-A*t) - B*y
应用积分因子exp(B*t)
得到等效值
(exp(B*t)*y(t))' = A*N0*exp((B-A)*t)
积分得到
exp(B*t)*y(t) = A*N0*exp((B-A)*t)/(B-A) + C
y(t) = A*N0*exp(-A*t)/(B-A) + C*exp(-B*t)
这正是求解器的计算结果。
您是否插入了 C1 = 0,并期望得到 y(0) = 0 的解?它不是这样工作的。 C1 是公式中的任意常数,将其设置为 0 并不能保证当 t = 0 时表达式的计算结果为 0。
这是一个正确的方法,循序渐进
sol1 = sp.dsolve(eq, daughter)
这个 returns 一个 Piecewise 因为不知道两个 lambda 是否相等:
Eq(N2(t), (C1 + N_0*lambda_1*Piecewise((t, Eq(lambda_2, lambda_1)), (-exp(lambda_2*t)/(lambda_1*exp(lambda_1*t) - lambda_2*exp(lambda_1*t)), True)))*exp(-lambda_2*t))
我们可以澄清他们不是:
sol2 = sol1.subs(Eq(lambda_2, lambda_1), False)
得到
Eq(N2(t), (C1 - N_0*lambda_1*exp(lambda_2*t)/(lambda_1*exp(lambda_1*t) - lambda_2*exp(lambda_1*t)))*exp(-lambda_2*t))
接下来,我们需要C1,使得当t = 0时右侧变为0。所以,取右侧,将0代入t,求解C1:
C1 = Symbol('C1')
C1val = solve(sol2.rhs.subs(t, 0), C1, dict=True)[0][C1]
(没有必要包含 dict=True
但我喜欢它,因为它强制统一输出 solve:它是一个字典数组。)
顺便说一下,C1val 现在是 N_0*lambda_1/(lambda_1 - lambda_2)
。放入:
sol3 = sol2.subs(C1, C1val).simplify()
你有它:
Eq(N2(t), N_0*lambda_1*(exp(lambda_1*t) - exp(lambda_2*t))*exp(-t*(lambda_1 + lambda_2))/(lambda_1 - lambda_2))
该表达式等同于 N_0*lambda_1*(exp(-lambda_2*t) - exp(-lambda_1*t))/(lambda_1 - lambda_2)
,尽管 SymPy 似乎不愿意在此处组合指数。
我正在使用 sympy.dsolve 求解衰减链的简单 ODE。 对于不同的衰减率(例如 lambda_1 > lambda_2),我得到的答案是错误的。代入 C1=0 后,我得到一个简单的指数
-N_0*lambda_1*exp(-lambda_1*t)/(lambda_1 - lambda_2)
而不是正确答案有:
(exp(-lambda_1*t)-exp(-lambda_2*t)).
我做错了什么? 这是我的代码
sp.var('lambda_1,lambda_2 t')
sp.var('N_0')
daughter = sp.Function('N2',Positive=True)(t)
stage1 = N_0*sp.exp(-lambda_1*t)
eq = sp.Eq(daughter.diff(t),stage1*lambda_1 - daughter*lambda_2)
sp.dsolve(eq,daughter)
您的微分方程是(使用较短的变量标识符)
y' = A*N0*exp(-A*t) - B*y
应用积分因子exp(B*t)
得到等效值
(exp(B*t)*y(t))' = A*N0*exp((B-A)*t)
积分得到
exp(B*t)*y(t) = A*N0*exp((B-A)*t)/(B-A) + C
y(t) = A*N0*exp(-A*t)/(B-A) + C*exp(-B*t)
这正是求解器的计算结果。
您是否插入了 C1 = 0,并期望得到 y(0) = 0 的解?它不是这样工作的。 C1 是公式中的任意常数,将其设置为 0 并不能保证当 t = 0 时表达式的计算结果为 0。
这是一个正确的方法,循序渐进
sol1 = sp.dsolve(eq, daughter)
这个 returns 一个 Piecewise 因为不知道两个 lambda 是否相等:
Eq(N2(t), (C1 + N_0*lambda_1*Piecewise((t, Eq(lambda_2, lambda_1)), (-exp(lambda_2*t)/(lambda_1*exp(lambda_1*t) - lambda_2*exp(lambda_1*t)), True)))*exp(-lambda_2*t))
我们可以澄清他们不是:
sol2 = sol1.subs(Eq(lambda_2, lambda_1), False)
得到
Eq(N2(t), (C1 - N_0*lambda_1*exp(lambda_2*t)/(lambda_1*exp(lambda_1*t) - lambda_2*exp(lambda_1*t)))*exp(-lambda_2*t))
接下来,我们需要C1,使得当t = 0时右侧变为0。所以,取右侧,将0代入t,求解C1:
C1 = Symbol('C1')
C1val = solve(sol2.rhs.subs(t, 0), C1, dict=True)[0][C1]
(没有必要包含 dict=True
但我喜欢它,因为它强制统一输出 solve:它是一个字典数组。)
顺便说一下,C1val 现在是 N_0*lambda_1/(lambda_1 - lambda_2)
。放入:
sol3 = sol2.subs(C1, C1val).simplify()
你有它:
Eq(N2(t), N_0*lambda_1*(exp(lambda_1*t) - exp(lambda_2*t))*exp(-t*(lambda_1 + lambda_2))/(lambda_1 - lambda_2))
该表达式等同于 N_0*lambda_1*(exp(-lambda_2*t) - exp(-lambda_1*t))/(lambda_1 - lambda_2)
,尽管 SymPy 似乎不愿意在此处组合指数。