Matlab中龙格-库塔(Runge-Kutta)方法原理及实现

2023-11-25 18:40

本文主要是介绍Matlab中龙格-库塔(Runge-Kutta)方法原理及实现,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

本文为转载ilovematlab中hyowinner大神的文章,向前辈致敬。

龙格-库塔(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);

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

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

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数)
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值,可以是向量,但是要注意维数
for ii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end

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

function dy=test_fun(x,y)
dy = zeros(3,1);%初始化列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) + y(3);
dy(3) = -0.51 * y(1) * y(2);
对该微分方程组用ode45和自编的龙格库塔函数进行比较,调用如下:
[T,F] = ode45(@test_fun,[0 15],[1 1 3]);
subplot(121)
plot(T,F)%Matlab自带的ode45函数效果
title('ode45函数效果')
[T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数,别忘记改变初始值的维数
subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果title('自编的   龙格库塔函数')


这篇关于Matlab中龙格-库塔(Runge-Kutta)方法原理及实现的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/424348

相关文章

pandas中位数填充空值的实现示例

《pandas中位数填充空值的实现示例》中位数填充是一种简单而有效的方法,用于填充数据集中缺失的值,本文就来介绍一下pandas中位数填充空值的实现,具有一定的参考价值,感兴趣的可以了解一下... 目录什么是中位数填充?为什么选择中位数填充?示例数据结果分析完整代码总结在数据分析和机器学习过程中,处理缺失数

Golang HashMap实现原理解析

《GolangHashMap实现原理解析》HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持高效的插入、查找和删除操作,:本文主要介绍GolangH... 目录HashMap是一种基于哈希表实现的键值对存储结构,它通过哈希函数将键映射到数组的索引位置,支持

Java学习手册之Filter和Listener使用方法

《Java学习手册之Filter和Listener使用方法》:本文主要介绍Java学习手册之Filter和Listener使用方法的相关资料,Filter是一种拦截器,可以在请求到达Servl... 目录一、Filter(过滤器)1. Filter 的工作原理2. Filter 的配置与使用二、Listen

Pandas使用AdaBoost进行分类的实现

《Pandas使用AdaBoost进行分类的实现》Pandas和AdaBoost分类算法,可以高效地进行数据预处理和分类任务,本文主要介绍了Pandas使用AdaBoost进行分类的实现,具有一定的参... 目录什么是 AdaBoost?使用 AdaBoost 的步骤安装必要的库步骤一:数据准备步骤二:模型

Pandas统计每行数据中的空值的方法示例

《Pandas统计每行数据中的空值的方法示例》处理缺失数据(NaN值)是一个非常常见的问题,本文主要介绍了Pandas统计每行数据中的空值的方法示例,具有一定的参考价值,感兴趣的可以了解一下... 目录什么是空值?为什么要统计空值?准备工作创建示例数据统计每行空值数量进一步分析www.chinasem.cn处

使用Pandas进行均值填充的实现

《使用Pandas进行均值填充的实现》缺失数据(NaN值)是一个常见的问题,我们可以通过多种方法来处理缺失数据,其中一种常用的方法是均值填充,本文主要介绍了使用Pandas进行均值填充的实现,感兴趣的... 目录什么是均值填充?为什么选择均值填充?均值填充的步骤实际代码示例总结在数据分析和处理过程中,缺失数

Java对象转换的实现方式汇总

《Java对象转换的实现方式汇总》:本文主要介绍Java对象转换的多种实现方式,本文通过实例代码给大家介绍的非常详细,对大家的学习或工作具有一定的参考借鉴价值,需要的朋友参考下吧... 目录Java对象转换的多种实现方式1. 手动映射(Manual Mapping)2. Builder模式3. 工具类辅助映

Go语言开发实现查询IP信息的MCP服务器

《Go语言开发实现查询IP信息的MCP服务器》随着MCP的快速普及和广泛应用,MCP服务器也层出不穷,本文将详细介绍如何在Go语言中使用go-mcp库来开发一个查询IP信息的MCP... 目录前言mcp-ip-geo 服务器目录结构说明查询 IP 信息功能实现工具实现工具管理查询单个 IP 信息工具的实现服

SpringBoot基于配置实现短信服务策略的动态切换

《SpringBoot基于配置实现短信服务策略的动态切换》这篇文章主要为大家详细介绍了SpringBoot在接入多个短信服务商(如阿里云、腾讯云、华为云)后,如何根据配置或环境切换使用不同的服务商,需... 目录目标功能示例配置(application.yml)配置类绑定短信发送策略接口示例:阿里云 & 腾

Windows 上如果忘记了 MySQL 密码 重置密码的两种方法

《Windows上如果忘记了MySQL密码重置密码的两种方法》:本文主要介绍Windows上如果忘记了MySQL密码重置密码的两种方法,本文通过两种方法结合实例代码给大家介绍的非常详细,感... 目录方法 1:以跳过权限验证模式启动 mysql 并重置密码方法 2:使用 my.ini 文件的临时配置在 Wi