从系列中:MATLAB中常微分方程的求解
克利夫·莫勒
ODE4采用了经典的龙格-库塔方法,龙格-库塔方法是过去100年来ODE4使用最广泛的数值方法。它的主要缺点是缺乏误差估计。在这里和后面的视频中会用到一个火焰增长的简单模型。
这是经典的龙格-库塔方法。在20世纪上半叶的手工计算和20世纪后半叶的数字计算机计算中,这是100多年来世界上最流行的数值方法。我怀疑它今天还在使用。
每步计算函数四次,第一次是在区间的开始。然后用它来步进到区间的中间,得到s2。然后用s2再次步进到区间的中间。
然后再次计算函数,得到s3。然后使用s3在整个时间间隔内逐步清除,得到s4。然后把这四个斜坡组合起来,把中间的两个更重,最后走一步。这是经典的龙格库塔方法。
这是我们的MATLAB实现。我们称之为ODE4,因为它每一步计算四次函数。相同的参数,向量y out。现在我们有四个斜坡——S1开始,S2中途,S3再次在中间,然后在右边S4。s1的1/6、s2的1/3、s3的1/3和s4的1/6为您提供了最后一步。这是经典的龙格库塔方法。
卡尔·伦格(Carl Runge)是一位相当杰出的德国数学家和物理学家,他于1895年与其他几位学者一起发表了这一方法。他发表了许多其他的数学论文,并且相当有名。
马丁·库塔独立发现了这种方法,并于1901年发表。他在其他方面几乎不出名。
我想讲一个简单的燃烧模型。因为模型有一些重要的数值性质。如果你点燃一根火柴,火球会迅速增长,直到达到临界大小。然后,由于球内部燃烧所消耗的氧气量与球表面的氧气量相平衡,剩余的氧气量就会达到这个尺寸。
这是无量纲模型。匹配的是一个球面,y是它的半径。y的立方项是球的体积。y的立方表示内部的燃烧。
球面与y的平方成正比。y的平方项表示了表面上的氧。关键参数,重要参数,是初始半径,y0。
半径从y0开始增长,直到y立方项和y平方项相互平衡,此时增长率为0。半径不再增长。我们在很长一段时间内融合在一起。我们积分的时间与初始半径成反比。这就是模型。
这是一个动画。我们从一个小火焰开始,一个小球形火焰。你会看到一个小半径。时间和半径显示在图的顶部。它开始增长。
当时间到50秒时,我们已经做了一半了。火焰就像爆炸一样,然后上升到半径,这时这两项相互平衡。火焰停止生长。它在这里还在增长,虽然你在这个规模上看不出来。
让我们为龙格-库塔做准备。微分方程是y ' = y²- y³。从0开始,临界初始半径是0.01。这意味着我们要积分到2 / y0到时间200。
我要选择500步的步长。我随便选一个。好的,现在我准备使用ODE4。我把结果存储在y中。
一直到1。我要把结果画出来。这就是我需要的t的值。这是情节。
现在,你可以看到火焰开始膨胀。它生长得相当缓慢。然后在时间间隔的中间,它爆炸了,并且快速上升,直到到达半径为1的地方,然后停留在这里。
这个过渡时期是相当狭窄的。我们将继续研究这个问题。这个过渡区域将会给数值方法带来挑战。
现在,我们刚刚经历了它。我们有一个步长h,这是我们任意选择的。我们刚刚生成了这些值。我们真的不知道这些数字有多准确。
它们看起来不错。但它们有多准确?这是关于经典龙格-库塔方法的关键问题。图中的值有多可靠?
我有四个练习供你参考。如果微分方程不包含y,那么这个解就是一个积分。龙格-库塔法成为一种经典的数值积分方法。如果您学习过这些方法,那么您应该能够识别这种方法。
号码。第二,求出y ' = 1 + y²,y(0) = 0的精确解。然后看看ODE4会发生什么,当你尝试在t从0到2的区间上解它。
第三点,如果间隔的长度不能被步长整除会怎样?例如,如果t final是pi,并且步长是0.1。不要试图解决这个问题。这只是固定步长的风险之一。
最后,练习四——研究初始半径为1/1000的火焰问题。对于t的哪个值,半径达到其最终值的90%?
你也可以从以下列表中选择一个网站:
选择中国站点(中文或英文)以获得最佳站点性能。其他MathWorks国家/地区网站未针对您所在地的访问进行优化。