二阶偏微分方程组 龙格库塔法_四阶龙格库塔法

(建议阅读原文)预备知识 中点法解常微分方程(组)龙格库塔法是一类数值解微分方程的算法, 其中较常见的是四阶龙格库塔法. 这里不进行推导, 仅仅给出公式如下(

的定义类比式 5 )

其中

由以上两式, 不难把该算法拓展到方程组的情况. 对于
元微分方程组

我们可以把该式记为矢量函数的形式

现在我们仅需要把式 1 和 式 2 中的所有
都变为
维列矢量
即可将微分方程拓展为微分方程组.
例程
   到目前为止, 我们每求一个微分方程的数值解都要重新写一次程序, 对于一些较为复杂的算法这样做效率较低. 我们这里不妨把四阶龙格库塔法写到一个单独的函数文件 odeRK4.m 中, 当我们要解某个特定的方程时, 只需把式 2 中的
作为自变量输入即可解出

代码 1:odeRK4.m

8d2d1fc766d7295a49fc48fb5e052465.png


   我们先来看第 1 行的函数声明, 输入变量中,f 是式 4 中

的函数句柄,
tspan 是一个 2×1 的列矢量, tspan(1) 是初始时间, tspan(2) 是终止时间, Y0 是一个列矢量, Y0(ii) 是第 ii 个因变量的初始值, Nt
的个数,
tspan 定义的时间区间被等分为 Nt - 1 个小区间. 因变量中, Y 的行数是因变量的个数, 列数是 Ntt 是一个行矢量, 由第 6 行定义, Y(ii, jj) 就是第 ii 个变量在 t(jj) 时刻的值. 第 5 行把初值 Y0 赋给 Y 的第 1 列, 第 8-14 行的循环根据式 1 和式 2 的矢量形式由 Y 的第 ii 列(
)求第
ii+1 列(
).

   我们先来用这个函数来计算 “天体运动的简单数值计算” 中的问题. 我们令因变量
的四个分量依次为一阶方程组(式 4 )

中的
. 程序代码如下

代码 2:keplerRK4.m

0e589671e0092778004320ab15db8ae2.png

091aa55f2387b16766f5d6388792835e.png

运行结果如图 1 所示.

08399112c2ecaf6808b503e0827c76da.png
图 1:运行结果


   我们先来看函数 fun (20 行), 这个函数就相当于式 5 . 第一个输入变量 Y 是一个列矢量, 是

的值, 第二个输入变量是
, 但由于式 5 中没有出现
, 我们用波浪线代替. 第三个输入变量是参数
, 即万有引力常数和中心天体质量之积. 输出变量
Y1 是一个列矢量, 是
的值.

   再来看主函数 KeplerRK4 (第 1 行), 参数设定中除了步数从 4000 变为了 100, 其他都和 “天体运动的简单数值计算” 中的程序一样, 然而这里运行结果却精确得多(曲线几乎闭合), 可见这种算法的优越性.
   主函数第 10 行中将 fun(Y, t, GM) 变为函数句柄 f(Y, t), 这样 GM 就可以在 “参数设定” 中设置, 而不用在 fun 函数内部设置. 第 11 行调用了上文中的 odeRK4 函数解方程组, 由于我们在画图时不需要用 t, 所以把第二个输出变量改为波浪线.
表情包
插入表情
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符
相关推荐
©️2020 CSDN 皮肤主题: 游动-白 设计师:白松林 返回首页