声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4249|回复: 3

[编程技巧] Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

[复制链接]
发表于 2017-1-19 11:51 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分,它在积分区间多预计算出几个点的斜率,然后进行加权平均,用做下一点的依据,从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法,如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作:y'=f(x,y),使用差分概念。

(Yn+1-Yn)/h= f(Xn,Yn)推出(近似等于,极限为Yn')
Yn+1=Yn+h*f(Xn,Yn)
另外根据微分中值定理,存在0<t<1,使得
Yn+1=Yn+h*f(Xn+th,Y(Xn+th))
这里K=f(Xn+th,Y(Xn+th))称为平均斜率,龙格库塔方法就是求得K的一种算法。

利用这样的原理,经过复杂的数学推导(过于繁琐省略),可以得出截断误差为O(h^5)的四阶龙格库塔公式:
K1=f(Xn,Yn);
K2=f(Xn+h/2,Yn+(h/2)*K1);
K3=f(Xn+h/2,Yn+(h/2)*K2);
K4=f(Xn+h,Yn+h*K3);
Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6);

所以,为了更好更准确地把握时间关系,应自己在理解龙格库塔原理的基础上,编写定步长的龙格库塔函数,经过学习其原理,已经完成了一维的龙格库塔函数。

仔细思考之后,发现其实如果是需要解多个微分方程组,可以想象成多个微分方程并行进行求解,时间,步长都是共同的,首先把预定的初始值给每个微分方程的第一步,然后每走一步,对多个微分方程共同求解。想通之后发现,整个过程其实很直观,只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数:

  1. function [x,y]=runge_kutta1(ufunc,y0,h,a,b)
  2. %参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)
  3. n=floor((b-a)/h);%求步数

  4. x(1)=a;%时间起点

  5. y(:,1)=y0;%赋初值,可以是向量,但是要注意维数

  6. for ii=1:n

  7. x(ii+1)=x(ii)+h;

  8. k1=ufunc(x(ii),y(:,ii));

  9. k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);

  10. k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);

  11. k4=ufunc(x(ii)+h,y(:,ii)+h*k3);

  12. y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
  13. %按照龙格库塔方法进行数值求解
  14. end
复制代码


调用的子函数以及其调用语句:

  1. function dy=test_fun(x,y)
  2. dy = zeros(3,1);%初始化列向量
  3. dy(1) = y(2) * y(3);
  4. dy(2) = -y(1) + y(3);
  5. dy(3) = -0.51 * y(1) * y(2);
复制代码


对该微分方程组用ode45和自编的龙格库塔函数进行比较,调用如下:
  1. [T,F] = ode45(@test_fun,[0 15],[1 1 3]);
  2. subplot(121)
  3. plot(T,F)%Matlab自带的ode45函数效果
  4. title('ode45函数效果')
  5. [T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数,别忘记改变初始值的维数
  6. subplot(122)
  7. plot(T1,F1)%自编的龙格库塔函数效果
  8. title('自编的   龙格库塔函数')
复制代码


运行结果如下:
203.png


转自:matlab中文网
作者:hyowinner
回复
分享到:

使用道具 举报

发表于 2017-1-19 12:52 | 显示全部楼层
好贴 多谢分享
发表于 2017-10-8 17:10 | 显示全部楼层
多谢楼主分享
发表于 2024-6-9 15:39 | 显示全部楼层
一直对其复杂的数学推导过程感兴趣
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-10 12:59 , Processed in 0.072902 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表