在MATLAB中求解ode, 8:方程组
从系列中:用MATLAB求解ode
含高阶导数的常微分方程重写为只含一阶导数的矢量系统。以经典的范德波尔非线性振荡器为例。随着参数的增加,VdP方程变得僵硬。
许多数学模型都涉及到高阶导数。但是MATLAB ODE求解器只适用于一阶常微分方程组。所以我们要重写模型只包含一阶导数。我们来看看如何用一个非常简单的模型,谐振子。
X ' ' + X = 0。这涉及到二阶导数。为了把它写成一阶方程组,我们引入了向量y,这是一个有两个分量的向量,x和x的导数。
我们只是改变了符号,让y有两个分量,x和x '。那么y的导数就是包含x和x ' '的向量。微分方程变成y ' + y = 0。看到我们是如何重写这个微分方程的了吗。y '是x ' '的导数?
一旦你做到了这一点,其他的事情就很容易了。向量系统现在是y1 y '是y2 - y1。第一个分量是y ' = y2。这就是说第一个分量的导数是第二个分量。这是微分方程本身。y ' = - y1是谐振子的微分方程。
当我们把它写成MATLAB的自主函数时,这就是自主函数。f是t和y的独立函数,它不依赖于t,首先它是一个向量,一个列向量。f的第一项是y2。第二个分量是-y。
第一个分量只是一个符号问题。所有的内容都在第二个分量中,它表示微分方程。
对于某些初始条件,假设初始条件是x (0) = 0 x '(0) = 1。用向量y表示,y的第一个分量是0。第二个分量是1。
或者用向量表示,初始向量是。这意味着它们的解是sint和cost,当我们在MATLAB中写出初始条件时,它是列向量。
让我们打开MATLAB命令窗口。这是微分方程。y ' = y2。y ' = -y。这是谐振子。
我们从0到2积分,因为它们是三角函数。我要求输出的步长为2 / 36,相当于每10度,就像机场的跑道一样。
我需要一个初始条件。Y0 0 = 0。我需要一个列向量,来表示这两个分量。我将使用ODE45,如果我不带输出参数调用它,微分方程的ODE45 (f, t)张成时间区间,y0是初始条件。
如果我不带输出参数调用ODE45,它就会自动绘制解决方案。这里我们得到了cost从1开始,sint从0开始的图像。
现在如果我回到命令窗口,并要求捕获t和y的输出,我就会得到输出的向量。37步,向量t,两个分量y,两列包含正弦和余弦。现在我可以把它们画在相平面上。
画出ya / y2。如果我对坐标轴进行正则化,就会得到一个漂亮的圆图,每10度就有一个小圆,就像机场的跑道一样。
范德堡尔振荡器于1927年由一位荷兰电气工程师引入,用于模拟真空管电路中的振荡。它有一个非线性阻尼项。从那以后,它被用来模拟各种领域的现象,包括地质学和神经学。
它表现出混沌的行为。我们对数值分析感兴趣是因为,随着参数mu的增大,问题变得越来越复杂。为了将它写成一阶系统,以便与MATLAB ODE求解器一起使用,我们引入向量y,包含x和x '。
所以y ' = x '和x ' '然后微分方程被写成y '的第一个分量是y2。然后微分方程写成y的第二分量。
这是非线性阻尼项- y1。当为0时,它就变成了谐振子。这是匿名函数。
让我们用范德堡尔振荡器做一些实验。首先,我要定义mu的值。取一个适中的值。等于100。现在有了集合,我可以定义匿名函数。
它包含一个mu的值。这是f的第二个分量中的范德堡尔方程,我将收集关于ODE求解器工作难度的统计数据。为此,我将使用ODE set,并告诉它我想打开stats。
我需要一个初始条件。现在我要运行ODE45来解决这个问题。我只指定一个起始值t,和一个最终值t ODE45会选择它自己的时间步长。我知道t的现价是460,它会对这两个周期的解进行积分。
现在我们可以看一段时间了。这需要很多步骤。随着步骤越来越多,速度开始变慢。现在它开始变得非常缓慢,因为它开始变硬。
我先把摄像机关掉,这样你就不用看这些步骤了。我们正试着升到460。快结束的时候我再把它打开。
好吧,摄像机已经坏了三分钟了。你可以看到我们已经取得了多大进展。我们离终点还远着呢。我要在这里按下停止按钮。然后回到命令窗口。
哦,我们还没讲到最后。让我退回到时间区间试试这个值。看看这是怎么回事。所以这仍然需要很多步骤。
但是我们可以,这大概是一个周期。我们会讲到最后。
我会让摄像机开着直到我们结束。这花了不到一分钟的时间。它走了将近15000步。让我们用硬解器来求解。
在那里。所以只用了半秒,只走了500步。这里有一个简单的刚度例子。让我们用刚性求解器来检验范德堡尔方程。让我们捕获输出并绘制它。
因为那个情节不是很有趣。我想用几种不同的方式画出来。再一次,我想回到——捕捉几个时期。
让我们画出一个电流分量作为时间的函数。就是这个。这是范德堡尔方程。
你可以看到它有一个初始的暂态,然后它进入周期性振荡,这里有非常陡峭的尖峰。即使是这个僵硬的求解者也在努力处理这些快速变化。这里有很多点,因为它选择了步长。
现在,让我们回到命令行,画相平面图。这是带阻尼的振子的相平面图。它不在圆附近,如果是0,它就在圆附近。
这就是范德堡尔振荡器的特性。
您也可以从以下列表中选择一个网站:
如何获得最佳的网站性能
选择中国站点(中文或英文)以获得最佳站点性能。其他MathWorks国家站点没有针对您所在位置的访问进行优化。