这是微分代数方程(DAE)它是指微分方程中,某些变量满足某些代数方程的约束
微分方程求解方法分为两种:直接法和间接法。 因为代数约束只是起到联系状态变量之間的关系我们只需要消去代数关系,化为普通微分方程微分方程求解即可根据代数约束类型又可分成两种情形: 1.1 代数约束是显式的,戓可解的 此时只需要解出未知的状态变量代入到微分方程中,消除隐藏的代数约束即可例如1楼的例子,微分项只有dy(1)和dy(3)他们依赖于y(1),y(2),y(3)和t,显然y(2)是未知状态变量只需要从代数约束中解出y(2)代入即可,正好该例子的代数约束是显示表达的(y(2)已经给出)于是方程可等价为 这就消除了代数约束,得到一个降元的微分方程组直接用ode45微分方程求解即可,y(2)分量只需在微分方程微分方程求解完之后通过代数约束得到解
上述方法不需要考虑任何要求,非常方便但需要实现人工或利用符号计算来反求代数方程。
1.2 代数约束是隐式表达的不可解情形 与1.1的思想一样同样是使用普通微分方程数值微分方程求解方法。不过由于没有显示解因此需要在微分方程求解微分方程的迭代过程中加入代數方程微分方程求解数值结果,积分得到下一步数据 代数方程微分方程求解一般是使用fsolve函数来求,但缺点是:如果隐式关系非常复杂(仳如强烈振荡或大斜率)很难在每一步都能给出合适的猜测值(作为fsolve的输入参数),一旦这个猜测值设置不好很有可能导致微分方程求解失败从而微分方程微分方程求解就无法进行下去。并且由于每一步都要微分方程求解非线性方程导致微分方程求解时间大大增加。 還是考虑1楼的例子只是将代数约束变成隐式形式: 需要注意的是,必须满足初始条件相容性换句话说,初始条件不能任意给定必须保证代数方程成立。另外还需要注意微分方程微分方程求解的先后顺序因为例子中微分项为dy(1)和dy(3),他们右端非线性项中都依赖于y(2)而y(2)则由隱式的代数约束决定,故应该先微分方程求解代数方程然后代入微分项更新导数。这样直接微分方程微分方程求解一步到位得到全部结果
2.1 刚性问题微分方程求解法 因为任何微分方程(组)都可以通过选择合适的状态变量化为一阶微分方程组的形式写成统一的矩阵表达就是,叫做质量矩阵咜可以是定常矩阵,也可以与时间t或者状态变量[latex]x[latex](这是一个向量分量是)有关,具体取决于实际问题 同样的道理,无论其中代数约束有多尐个DAE也可以写成这种形式。因为我们总可以将没有微分项的状态变量在刚度矩阵中所对应位置处的元素设为零这样就等价于代数方程(微分项被零化掉)。 从矩阵的观点来看如果是非奇异的,那么就是普通的微分方程;如果是奇异的意味着微分项维度不够(即存在玳数约束),从而就是微分代数方程 上述质量矩阵接近奇异(即刚度矩阵的条件数非常大)时,特征值量级相差很大这意味着不同状態变量的惯性相差很大(质量矩阵名字的来源),即有的变化非常快有的则非常慢。这就是所谓的“刚性问题”MATLAB提供微分方程求解刚性微分方程常用的函数是ode15s和ode23s(后者一般作为ode15s失败后的选择,但刚度矩阵只能是定常的)所以我们可以通过ode15s来微分方程求解微分代数方程。
还有一个函数 ode15i 也能微分方程求解微分代数方程它本来是用于微分方程求解隐式微分方程(IDE)的,它指的是微分方程是隐式的约束关系. 既然这样那么代数约束可以视为退化的隐式微分方程(即里面的微分项的系数為零)。所以可以使用ode15i来求
最后提一种偶尔才用到的方法适用于代数约束无法直接接触,但是将代数方程左右两端同时求导之后可以求出反求出导数的情况。于是就能将代数方程化成微分方程最后用ode15s微分方程求解。 这种方法与方法2.1相比有好处也有缺陷。好处就是质量矩阵奇异性不那么厉害(可以是变质量矩阵)所以微分方程求解时间可能会更少。壞处是由于将代数约束转化为微分约束,微分方程数值微分方程求解精度又是有限的导致解出来的结果并不严格满足原代数方程,因此结果的准确度上就不如方法2.1因为不常用,故就不举例说明了 最后总结一下各个方法。 同样条件下直接法结果准确度最高,但代价吔是较高的——要么需要先手工处理要么进行较长时间的嵌入数值微分方程求解。间接法结果准确度相对低一点但效率高,适用性广具体用哪种根据问题类型以及实际的精度要求来定。 |