我对 Mathematica NDSolve 中的约束有疑问
如果我有两个线性微分方程
NDSolve[x'[t]= A x[t]+B y[t],
y'[t]= C x[t]+D y[t],
x[0]==0.002, y[0]==0.005, {x,y}, {t,0.10000}]
我怎样才能施加约束
x[t]+y[t] == 1
对于任何时间 t
谢谢
我对 Mathematica NDSolve 中的约束有疑问
如果我有两个线性微分方程
NDSolve[x'[t]= A x[t]+B y[t],
y'[t]= C x[t]+D y[t],
x[0]==0.002, y[0]==0.005, {x,y}, {t,0.10000}]
我怎样才能施加约束
x[t]+y[t] == 1
对于任何时间 t
谢谢
答案是一般情况下你不能。这不是 Mathematica 的限制,而是数学的限制。
您正在尝试同时求解三个方程:
x' = Ax+By
y' = Cx+Dy
1 = x + y
我将使用小写字母a,b,c,d
来避免 Mathematica 出现问题。
正如贝利撒留建议的那样,您不需要手动进行替换。您可以简单地选择其中任何两个并要求 Mathematica 解决:
a = 1; b = 2; c = 3; d = 4; x0 = 0.2; y0 = 0.8;
NDSolve[{x'[t] == a x[t] + b y[t], y'[t] == c x[t] + d y[t],
x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, 0, 1000}];
Plot[Evaluate[ReplaceAll[{x[t], y[t]}, %]], {t, 0, 0.01}]
NDSolve[{x'[t] == a x[t] + b y[t], x[t]+y[t]==1,
x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, 0, 1000}];
Plot[Evaluate[ReplaceAll[{x[t], y[t]}, %]], {t, 0, 0.01}]
NDSolve[{y'[t] == c x[t] + d y[t], x[t]+y[t]==1,
x[0] == x0, y[0] == y0}, {x[t], y[t]}, {t, 0, 1000}];
Plot[Evaluate[ReplaceAll[{x[t], y[t]}, %]], {t, 0, 0.01}]
请注意,要使其中任何一个起作用,您将需要指定常量a,b,c,d,x0,y0
(并且您需要为第二个x0
和y0
第三个求和为 1)。
正如你所问的,我正在使用 NDSolve,尽管其中任何一个也可以通过DSolve
(有或没有常量的固定值)来解决:
Clear[a,b,c,d,x0,y0];
DSolve[{x'[t] == a x[t] + b y[t], y'[t] == c x[t] + d y[t],
x[0] == x0, y[0] == y0}, {x[t], y[t]}, t]
DSolve[{x'[t] == a x[t] + b y[t], x[t]+y[t]==1,
x[0] == x0, y[0] == 1-x0}, {x[t], y[t]}, t]
DSolve[{y'[t] == c x[t] + d y[t], x[t]+y[t]==1,
x[0] == x0, y[0] == 1-x0}, {x[t], y[t]}, t]
但是,非常重要的是要注意,对于三个函数方程和两个函数,您不能同时求解所有这些方程。这就是为什么生成的三张图NDSolve
和生成的三组解DSolve
会有所不同。请注意,在后两种情况下,我必须通过给出1-x0
初始值来明确强制 x[0]+y[0]=1。
因为你有三个方程,你要么有一个不一致的系统,要么是超定的。如果您尝试同时解决所有三个问题,Mathematica 将给出错误:
NDSolve::overdet:
"There are fewer dependent variables, {x[t],y[t]}, than equations, so the
system is overdetermined."
(类似的DSolve
。)
我的选择a=1,b=2,c=3,d=4
使它不一致(因此我的三个图表不同意)。你可以弄乱Manipulate
一些东西来看看它们是如何关联的,但是我有两种基本的方法来确定如何(如果有的话)强制解决方案对所有三个方程保持一致:
只是区分1=x+y
。这给了你0=x'+y'
这意味着你有整个限制ax+by=cx+dy
,换句话说y=(c-a)/(b-d)x
。这实际上与您原来的约束不一致x+y=1
。
考虑 DEsx'=ax+by
和的线性系统y'=cx+dy
。为了使解成为一条直线,它必须在特征向量的方向上并且必须通过原点。您的约束x+y=1
是一条不符合这些标准的线,因此两个 DE 的解都不会满足您的代数方程。
PS我认为您的原始代码中有错字。你可能想要{t,0,10000}
是的?
尝试例如:
DSolve[{x'[t] == A x[t] + B (1 - x[t]), x[0] == 2/10}, x[t], t]