MATLAB求解ode, 3:经典龙格库塔,ODE4
从系列中:用MATLAB求解ode
ODE4实现了经典的龙格-库塔方法,这是过去100年来ODE4最广泛使用的数值方法。它的主要缺点是缺乏误差估计。火焰生长的一个简单模型就是一个例子,在这里和后面的视频中都会用到。
这是经典的龙格-库塔方法。到目前为止,这是世界上最流行的数值方法,在20世纪上半叶的100多年里用于手工计算,然后在20世纪下半叶用于数字计算机的计算。我怀疑它现在还在使用。
每个步骤计算函数四次,第一次是在区间的开始。然后用这个步进到区间的中间,得到s2。然后使用s2再次进入间隔的中间。
然后对函数求值得到s3。然后用s3在区间内进行步进清理,得到s4。然后把这四个斜坡组合起来,中间的两个更重,来做最后一步。这就是经典的龙格-库塔方法。
这是我们的MATLAB实现。我们称它为ODE4,因为它每步计算4次函数。同样的参数,向量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的区间上解它。
第三,如果间隔的长度不能被步长整除会怎样?例如,如果tfinal是pi,步长是0.1。不要试图解决这个问题。这只是固定步长的危害之一。
最后,练习四,研究初始半径为千分之一的火焰问题。当t的值是多少时半径会达到最终值的90% ?
您也可以从以下列表中选择一个网站:
如何获得最佳的网站性能
选择中国站点(中文或英文)以获得最佳站点性能。其他MathWorks国家站点没有针对您所在位置的访问进行优化。