在MATLAB中求解ode, 3:经典龙格-库塔,ODE4
从系列:用MATLAB求解ode
ODE4实现了经典的龙格-库塔方法,这是在过去100年里最广泛使用的ODE4数值方法。它的主要缺点是缺乏误差估计。火焰生长的一个简单模型在这里和后面的视频中都有使用。
这是经典的龙格-库塔法。在20世纪上半叶的100多年里,这是世界上最流行的手工计算方法,在20世纪下半叶,这是在数字计算机上最流行的计算方法。我怀疑它现在还在使用。
每一步计算函数四次,第一次在区间的开始处。然后用它进入区间的中间,得到s2。然后再用s2步进到区间的中间。
然后对这个函数求值得到s3。然后使用s3清除间隔,得到s4。然后把这四个斜坡组合起来,把中间的两个斜坡加重,完成最后一步。这就是经典的龙格-库塔法。
这是我们的MATLAB实现。我们称它为ODE4,因为它每步计算四次函数。同样的,向量y出。现在我们有四个斜率,s1在开始,s2在中间,s3在中间,然后s4在右边。1/6 s1 1/3 s2 1/3 s3 1/6 s4就是最后一步。这就是经典的龙格-库塔法。
卡尔·朗格是一位相当杰出的德国数学家和物理学家,他在1895年和其他几个人一起发表了这个方法。他还发表了许多其他的数学论文,相当有名。
马丁·库塔独立发现了这种方法,并于1901年发表。他在其他任何方面都不太出名。
我想研究一个简单的燃烧模型。因为这个模型有一些重要的数值性质。如果你点燃一根火柴,火球迅速增长,直到达到临界大小。然后残骸保持这个尺寸,因为在球的内部燃烧消耗的氧气量平衡了通过表面的氧气量。
这是无量纲模型。火柴是一个球体,y是它的半径。y³项是球体的体积。y³表示内部的燃烧。
球面与y的平方成正比。y²项表示通过表面的氧气。关键参数,重要参数,是初始半径y0 y0。
半径从y0开始增长,直到y立方和y平方相互平衡,在这一点上增长率为0。半径不再增长。我们经过了很长时间的整合。积分时间与初始半径成反比。这就是模型。
这是一个动画。我们从一个小的火焰开始,一个小的球形火焰。你会看到一个小半径。时间和半径显示在图的顶部。它开始生长了。
当时间到50时,我们已经完成了一半。火焰会爆炸,然后上升到半径为1的地方,这时两项相互平衡。火焰停止生长。在这里它还在轻微增长,尽管你看不到它的规模。
让我们为龙格-库塔设置这个。微分方程是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是,步长是0.1。不要试图解决这个问题。这只是固定步长的危害之一。
最后,练习四,研究初始半径为1/ 1000的火焰问题。对于t的什么值,半径达到最终值的90% ?
您也可以从以下列表中选择网站:
如何获得最佳的网站性能
选择中国网站(中文或英文)以获得最佳的网站表现。其他MathWorks国家网站没有针对从您的位置访问进行优化。