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