【算法基础学习 6】龙格库塔法 求微分方程

龙格库塔法


龙格库塔法对于一个工科学生来说估计并不陌生,因为我们做项目的时候或多或少会接触到微分方程,解微分方程的数值解是就需要用到龙格库塔方程了。我第一次接触到该方法是对四轴姿态角进行四元数解算的时候,当时用的是一阶龙格库塔,当时的理解就是这不就是泰勒展开式取一阶嘛,然后觉得二阶龙格库塔应该就是取二阶了,依次类推,所以当时没有深入的了解,现在因为遇到这方面的问题所以重新把他学了一遍。


准确的说,高阶龙格库塔法是借用了泰勒展开式的思想,但是和泰勒展开式还是有区别的。下面我对龙格库塔法讲讲我的理解。

首先,龙格库塔的一般形式如下:

首先,我们可以知道一阶龙格库塔法的表达式

上述式子可以改写成如下式子

这里我们可以看成k是在步长为h情况下,函数的增长量,f(x,y)为函数y的斜率。这个式子其实比较好理解的。

然后我们来看二阶龙格库塔的推导过程:
二阶情况下一般表达式中就只有k1和k2了,二阶龙格库塔我们可以看成是一阶表达式的一个变形,由原来的一个k变成了k1,k2,用k1,k2来表示k,c1和c2相当于k1和k2所占的权重一样,这样理解是不是好一点。那k1和k2分别表示什么??其实k1表示xi点处的增长量,k2表示的是xi+a2*h处的增长量。在这基础上,因为c1、c2、a2和是未知量。我们需要将其求出来,下面是我的推导过程:


这样二阶龙格库塔的迭代式子就可以表达出来了,写完收工。
然而如果你想你的数值解更加精确的话,你就需要更高 阶次的龙格库塔迭代式,一般而言,四阶是最经典的,用的也是最多的。本来我想把四阶的参数也算一遍,但是工程量巨大,我放弃了。我尝试了一下3阶的参数推导,推导一半,我还是放弃了,我把我推到一半的过程也写出来,希望有人帮我弄出来,本人实在没有勇气。


这里把k1,k1k2,k2平方带入式子中,我感觉非常人能及的,高中学的解析几何感觉都不算什么??
这里有一个需要注意的地方,图片中我用红笔标注了,由于三阶龙格库塔的阶段误差为,意思就是多项式中存在h的次数大于等于4的都可以忽略,这可能会大大减少我们的计算量,但是还是好难。
最后给出一个四阶龙格库塔算迭代表达式:


 

  • 1
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
Matlab中使用龙格-库塔(Runge-Kutta)方法微分方程是一个常见的数值计算任务。以下是一个简单的4阶龙格-库塔(RK4)算法的示例代码,用于解一阶线性常微分方程: ```matlab function [t, y] = rk4(f, t0, y0, tfinal, h) % f: 微分方程的函数 handle (f(t, y) -> dy/dt) % t0: 初始时间 % y0: 初始值 % tfinal: 结束时间 % h: 时间步长 % 初始化 nsteps = floor((tfinal - t0) / h); t = zeros(1, nsteps + 1); % 存储时间点 y = zeros(size(y0)); % 存储每个时间步的解 % 龙格-库塔方法的核心步骤 for i = 1:nsteps k1 = h * f(t(i), y); k2 = h * f(t(i) + h/2, y + k1/2); k3 = h * f(t(i) + h/2, y + k2/2); k4 = h * f(t(i) + h, y + k3); % 更新解 y = y + (k1 + 2*k2 + 2*k3 + k4)/6; % 记录时间 t(i+1) = t(i) + h; end t(nsteps+1) = tfinal; % 设置最后一个时间点 ``` 这个函数`rk4`接受一个描述微分方程导数的匿名函数`f`作为输入,你需要自定义这个函数来表达你的具体方程。例如,如果你有一个方程`dy/dt = f(t, y)`,那么`f`会像这样定义:`f = @(t, y) ...;`。 你可以调用这个函数并传入你的方程、初始条件和参数,如下所示: ```matlab % 假设你有一个二阶线性方程 dy'' = -k*y k = 1; % 权重系数 f = @(t, y) [0; -k*y(1)]; % 示例使用 tspan = [0 10]; % 时间范围 y0 = [1; 0]; % 初始状态 [t, y] = rk4(f, tspan(1), y0, tspan(2), 0.1); ``` 运行后,`t`数组将包含时间点,`y`数组则包含了对应时间下的解。对于更复杂的微分方程或特定的需,可能还需要根据实际情况进行适当修改或扩展。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值