龙格-库塔法(Matlab实现)

news2024/11/26 22:31:29

四阶龙格-库塔法介绍

在各种龙格-库塔法当中有一个方法十分常用,以至于经常被称为“RK4”或者就是“龙格-库塔法”。该方法主要是在已知方程导数和初始值时,利用计算机的仿真应用,省去求解微分方程的复杂过程。

令初值问题表述如下:

则,对于该问题的RK4由如下方程给出:

其中:

这样,下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积所决定。该斜率是以下斜率的加权平均:

  • k1是时间段开始时的斜率;
  • k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn + h/2的值;
  • k3也是中点的斜率,但是这次采用斜率k2决定y值;
  • k4是时间段终点的斜率,其y值用k3决定。

当四个斜率取平均时,中点的斜率有更大的权值:

RK4法是四阶方法,也就是说每步的误差是h5阶,而总积累误差为h4阶。

注意上述公式对于标量或者向量函数(y可以是向量)都适用。

举例分析

以2022年国赛A题 波浪能最大输出功率设计 第一问为例(有些公式出现的很突兀,建议看一下优秀论文,更易理解)参考资料:2022高教社杯全国大学生数学建模竞赛A题论文展示(A001)

第一问建立的数学模型如下:

需要求解的是浮子与振子的位移和速度,标红表示已知参数

Matlab编程实现

第一题第(1)问完整代码如下:

% 第一题第(1)问
tic

clear;
clc;
close all;

%% 数据初始化
w = 1.4005; %角速度 入射波浪频率
T = (2*pi)/w*40; % 前40个波浪周期
h = 1e-4; % 步长
t = 0:h:T; % 生成自变量t的向量

mf = 4866; % 浮子的质量kg
mz = 2433; % 振子的质量kg
mfu = 1335.535; % 垂荡附加质量
S = mf + mfu;

C_x = 656.3616; % 兴波阻尼
rho = 1025; % 海水密度
g = 9.8; % 重力加速度
A = pi; % 圆柱的横截面积
zf = 6250; % 垂荡激励力振幅(N)
K_tang = 80000;% 弹簧刚度
B = rho*g*A;
V = (0.8*pi)/3;% 圆锥的体积

K_zn = 10000;% 阻尼系数
% 创建计算结果x,y,z,m的数组
N = length(t);
f_x = zeros(1,N); % 浮子的位移
z_x = zeros(1,N); % 振子的位移
f_v = zeros(1,N); % 浮子的速度
z_v = zeros(1,N); % 振子的速度

%% 四阶龙格库塔迭代求解
% k1 = f(tn,yn)
% k2 = f(tn+h/2,yn+h/2*k1)
% k3 = f(tn+h/2,yn+h/2*k2)
% k4 = f(tn+h/2,yn+h*k3)
% y(n+1) = yn + h/6*(k1+2*k2+2*k3+k4)
for i = 2:N
    tn = t(i-1);
    fx = f_x(i-1);
    zx = z_x(i-1);
    fv = f_v(i-1);
    zv = z_v(i-1);

    dfx1 = fv; % yn
    dzx1 = zv;
    dfv1 = (zf*cos(w*tn) + K_zn*(zv-fv) + K_tang*(zx-fx) - C_x*fv - B*fx)/S;
    dzv1 = -(K_zn*(zv-fv) + K_tang*(zx-fx))/mz;

    dfx2 = fv + h/2*dfv1; % yn+k1*h/2 (即v0 + a*t)
    dzx2 = zv + h/2*dzv1;
    dfv2 = (zf*cos(w*(tn+h/2)) + K_zn*(zv+h/2*dzv1-fv-h/2*dfv1) + ...
        K_tang*(zx+h/2*dzx1-fx-h/2*dfx1) - C_x*(fv+h/2*dfv1) - B*(fx+h/2*dfx1))/S;
    dzv2 = -(K_zn*(zv+h/2*dzv1-fv-h/2*dfv1) + K_tang*(zx+h/2*dzx1-fx-h/2*dfx1))/mz;

    dfx3 = fv + h/2*dfv2; % yn+h/2*k2
    dzx3 = zv + h/2*dzv2;
    dfv3 = (zf*cos(w*(tn+h/2)) + K_zn*(zv+h/2*dzv2-fv-h/2*dfv2) + ...
        K_tang*(zx+h/2*dzx2-fx-h/2*dfx2) - C_x*(fv+h/2*dfv2) - B*(fx+h/2*dfx2))/S;
    dzv3 = -(K_zn*(zv+h/2*dzv2-fv-h/2*dfv2) + K_tang*(zx+h/2*dzx2-fx-h/2*dfx2))/mz;

    dfx4 = fv + h*dfv3; % yn+h/2*k3
    dzx4 = zv + h*dzv3;
    dfv4 = (zf*cos(w*(tn+h)) + K_zn*(zv+h*dzv3-fv-h*dfv3) + ...
        K_tang*(zx+h*dzx3-fx-h*dfx3) - C_x*(fv+h*dfv3) - B*(fx+h*dfx3))/S;
    dzv4 = -(K_zn*(zv+h*dzv3-fv-h*dfv3) + K_tang*(zx+h*dzx3-fx-h*dfx3))/mz;

    f_x(i) = fx + h/6*(dfx1+2*dfx2+2*dfx3+dfx4);
    z_x(i) = zx + h/6*(dzx1+2*dzx2+2*dzx3+dzx4);
    f_v(i) = fv + h/6*(dfv1+2*dfv2+2*dfv3+dfv4);
    z_v(i) = zv + h/6*(dzv1+2*dzv2+2*dzv3+dzv4);
end

%% 画图
figure();
subplot('121');
plot(t,f_x,'r');
hold on;
plot(t,z_x-f_x,'b');
plot([0,T],[0,0]);
legend('浮子','振子');
xlabel('时间 s');
ylabel('位移 m');

subplot('122');
plot(t,f_v,'r');
hold on 
plot(t,z_v-f_v,'b');
plot([0,T],[0,0]);
legend('浮子','振子');
xlabel('时间 s');
ylabel('速度 m/s');

resfx=[];
resfv=[];
reszx=[];
reszv=[];
cout=0;
%遍历时间间隔为0.2,前40个波浪周期的垂荡位移和速度
start=find(t==0);
over=find(t==0.2);
dtt=over-start;
index=1;
for xx=0:0.2:T
    cout=cout+1;
    %查找对应时间所在的位置
    %记录位移速度
    resfx(cout)=f_x(index);
    reszx(cout)=z_x(index);
    resfv(cout)=f_v(index);
    reszv(cout)=z_v(index);
    index=index+dtt;
end
time =0:0.2:T;
result = [time',resfx',resfv',reszx',reszv'];
xlswrite('myresult-test',result);

toc

运行结果:

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/2063063.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

干货分享|如何使用SD插件进行老旧照片修复上色?

每个家庭都保存着一些温馨记忆的老照片。修复并给老照片上色曾经是一项难度颇大的技术活,现在有了AI技术的加持,使用Photoshop和SD插件,几分钟内就能让那些泛黄老旧的照片焕然一新。 打开一张老照片后按CtrlA快捷键以选取整个画布&#xff…

Linux命令抽象

linux命令都含有一定格式,有具体的语法。我们应用命令时一般需要按语法应用。 有些特殊命令,不遵从通用格式,应用时要格外注意。 命令很多,不需要都记住,但对命令有一个整体的了解,对快速使用命令、找到需…

MySQL运维学习(2):主从复制

1.什么是主从复制 *主从复制是指将主数据库的DDL和 DML操作通过二进制日志传到从库服务器中,然后在从库上对这些日志重新执行(也叫重做),从而使得从库和主库的数据保持同步。 MySQL支持一台主库同时向多台从库进行复制,从库也可以作为其他从…

Django 集成与扩展:Database Task Queue

文章目录 Django 连接 PostgreSQL安装 PostgreSQL 驱动配置 PostgreSQL 数据库更新 Django 设置确保证书文件的权限测试连接示例:完整的 settings.py 配置注意事项 Django 连接 MySQL安装 MySQL 驱动配置 MySQL 数据库更新 Django 设置运行迁移调试连接问题 Django …

独立站卖家投放Facebook广告的8个建议

在外贸领域,独立站结合Facebook广告投放是一种推动销售增长的关键策略,而结合自动化广告投放工具能使广告投放变得更加高效。以下是一系列针对外贸独立站卖家的Facebook广告投放建议,希望能帮助卖家更有效地利用这一平台,实现营销…

ISO 26262中的失效率计算:IEC 61709-Clause 17_Switches and push-buttons

概要 IEC 61709是国际电工委员会(IEC)制定的一个标准,即“电子元器件 可靠性 失效率的基准条件和失效率转换的应力模型”。主要涉及电学元器件的可靠性,包括失效率的基准条件和失效率转换的应力模型。本文介绍IEC 61709第十七章&…

四川财谷通,信息科技引领者!

在数字化浪潮席卷全球的今天,电子商务作为新经济形态的重要代表,正以前所未有的速度改变着我们的生活方式和消费习惯。四川财谷通信息技术有限公司,作为这一领域的佼佼者,凭借其深厚的技术底蕴与创新思维,在抖音小店这…

机房动环监控系统的主要功能@卓振思众

机房动环监控系统(Data Center Environmental and Monitoring System)是一种用于监测和管理数据中心或机房内部环境和设备状态的系统。其主要目的是确保机房设备在最佳环境条件下运行,从而提高系统的稳定性和安全性。以下是【卓振思众】机房动…

QLibrary的load失败(0x000000c1)

前言 用vs加载dll库是没有问题&#xff0c;移植到qt creator开发却加载失败。 #include <QLibrary>void LoadDll() {QString appPath QCoreApplication::applicationDirPath();QString strLibFile appPath "/Pay.dll";QLibrary *m_pLib nullptr;if (QFile:…

行业标杆 | 澳鹏Appen入选“2024年中国AI基础数据服务研究报告”

AI基础数据服务可加速高质量数据的获取与标注&#xff0c;推动AI算法的创新与持续优化&#xff0c;是AI产业发展的重要支撑。艾瑞咨询近日发布《2024年中国AI基础数据服务研究报告》&#xff0c;深度剖析了当前AI数据行业的挑战和机遇&#xff0c;并前瞻预测了未来趋势。作为AI…

2.3.2存储修改调整

如果使用的是云存储&#xff0c;错误提示&#xff1a;这个点击生成海报&#xff0c;直接提示 二维码生成失败 修改方法路径&#xff1a;crmeb\services\QrcodeService.php 增加代码&#xff1a; (string) 2. 本地存储修改 &#xff1a; //return $this->setError(‘请检…

com.alibaba.fastjson.JSONException: unclosed string : 

场景: 解析json字符串到java对象中报错 FinanceDownLoadFileDto financeDownLoadFileDto JSON.parseObject(line, FinanceDownLoadFileDto.class); 分析: 这不用想,一定是json格式问题 ,但是我去核对了几次文本中的json格式是正确的,因为我是复制粘贴到代码中的,只有考虑是…

“失业程序员跑滴滴求生,意外踏入AI绘画新天地:一个家庭的逆境转机故事“

我叫李明泽&#xff0c;一名在IT行业摸爬滚打多年的程序员。在这个看似光鲜的职业背后&#xff0c;却隐藏着无数的心酸与无奈。曾几何时&#xff0c;我以为我会一直在这个行业稳稳当当&#xff0c;但现实却给了我一记响亮的耳光。 一、就业市场的寒冬 随着互联网行业增速放缓&a…

电脑回收站数据怎么恢复回来 回收站怎么恢复半年前的文件

回收站是电脑一项非常重要的功能。有些小伙伴在操作电脑的时候&#xff0c;可能会不小心将一些重要的文件资料误删除&#xff0c;这些误删除的文件资料&#xff0c;不会彻底的被删除&#xff0c;而是会暂时存储在回收站中&#xff0c;在一定程度上可以保证文件资料的“安全”。…

【机器学习-监督学习】双线性模型

【作者主页】Francek Chen 【专栏介绍】 ⌈ ⌈ ⌈Python机器学习 ⌋ ⌋ ⌋ 机器学习是一门人工智能的分支学科&#xff0c;通过算法和模型让计算机从数据中学习&#xff0c;进行模型训练和优化&#xff0c;做出预测、分类和决策支持。Python成为机器学习的首选语言&#xff0c;…

UDS 诊断 - ReadScalingDataByIdentifier(按标识符读取换算数据)(0x24)服务

UDS 诊断服务系列文章目录 诊断和通信管理功能单元 UDS 诊断 - DiagnosticSessionControl&#xff08;诊断会话控制&#xff09;&#xff08;0x10&#xff09;服务 UDS 诊断 - ECUReset&#xff08;ECU重置&#xff09;&#xff08;0x11&#xff09;服务 UDS 诊断 - SecurityA…

楼宇如何打造一个高效的智慧公厕系统

在现代化的楼宇中&#xff0c;公共设施的完善程度直接关系到人们的使用体验和楼宇的整体品质。其中&#xff0c;公厕作为必不可少的公共设施之一&#xff0c;其重要性不容忽视。打造一个高效的智慧公厕系统&#xff0c;不仅能够为使用者提供更加舒适、便捷的如厕环境&#xff0…

误删?损坏?SD卡数据恢复全攻略,让你的数据起死回生!

现在这年头&#xff0c;SD卡就像是我们数字生活的小助手。不管是拍照记录生活的人&#xff0c;还是玩无人机的高手&#xff0c;或者是上班经常传文件的白领&#xff0c;SD卡里都存着我们的重要信息。但是&#xff0c;万一这些信息出点问题&#xff0c;比如不小心删了、文件坏了…

云手机解决了TikTok哪些账号运营难题?

随着社交媒体的蓬勃发展&#xff0c;TikTok作为一款风靡全球的短视频应用&#xff0c;成为许多个人和企业进行品牌推广、内容创作的首选平台。然而&#xff0c;随之而来的是TikTok账号运营的一系列难题。本文将深入探讨云手机是如何解决这些难题的。 1、多账号运营的便捷性&…

怎么选择数据恢复精灵,如何恢复数据?

作为一名财务人员&#xff0c;每天打交道的除了账本就是各种电子表格和财务报告。数据的重要性不言而喻&#xff0c;一旦发生意外丢失&#xff0c;那可真是让人焦头烂额。今天给大家分享下三款优质数据恢复工具。 首先&#xff0c;我们得知道&#xff0c;数据丢失的原因多种多样…