四阶龙格库塔法求解二元二阶常微分方程

ops/2025/1/21 20:47:57/

龙格库塔法(Runge-Kutta methods)是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。在工程领域应用广泛,可用于求解复杂系统的运动方程等问题。
这里采用matlab程序编写代码实现龙格库塔法对于二元二阶常微分方程的求解。

{ x ¨ + x ˙ + x + 2 y ¨ = − c o s t , x ¨ + y ¨ = − s i n t − c o s t , x ( 0 ) = 0 , x ˙ ( 0 ) = 1 , y ( 0 ) = 1 , y ˙ ( 0 ) = 0 , \left\{ \begin{array}{lr} \ddot{x}+\dot{x}+x+2\ddot{y}=-\rm cos\it t, & \\ \ddot{x}+\ddot{y}=-\rm sin\it t-\rm cos\it t, \\ x(0)=0, \\ \dot{x}(0)=1, \\ y(0)=1, \\ \dot{y}(0)=0, \\ \end{array} \right. x¨+x˙+x+2y¨=cost,x¨+y¨=sintcost,x(0)=0,x˙(0)=1,y(0)=1,y˙(0)=0,
求当 t = ( 0 , 40 ) t=(0,40) t=(0,40)的时候 x x x y y y的值。
该方程有精确解
{ x = s i n t y = c o s t \left\{ \begin{array}{lr} x=\rm sin\it t \\ y=\rm cos\it t \\ \end{array} \right. {x=sinty=cost
下面我们通过四阶龙格库塔法进行求解:
首先设 { u 1 = x , u 2 = x ˙ w 1 = y , w 2 = y ˙ \left\{ \begin{array}{lr} u1=x,u2=\dot{x}\\ w1=y,w2=\dot{y} \end{array} \right. {u1=x,u2=x˙w1=y,w2=y˙
则上式变换为:
{ f 1 ( t ) = u 1 ˙ = u 2 f 2 ( t ) = u 2 ˙ = x ¨ = x ˙ + x − 2 s i n t − c o s t f 3 ( t ) = w 1 ˙ = w 2 f 4 ( t ) = w 2 ˙ = y ¨ = − x ˙ − x + s i n t \left\{ \begin{array}{lr} f1(t)=\dot{u1}=u2 \\ f2(t)=\dot{u2}=\ddot{x}=\dot{x}+x-2\rm sin\it t-\rm cos\it t \\ f3(t)=\dot{w1}=w2\\ f4(t)=\dot{w2}=\ddot{y}=-\dot{x}-x+\rm sin\it t \end{array} \right. f1(t)=u1˙=u2f2(t)=u2˙=x¨=x˙+x2sintcostf3(t)=w1˙=w2f4(t)=w2˙=y¨=x˙x+sint
则有四个子函数

function output = f1(x,u1,u2,w1,w2)
output = u2;
function output = f2(x,u1,u2,w1,w2)
output = u2+u1-2*sin(x)-cos(x);
function output = f3(x,u1,u2,w1,w2)
output = w2;
function output = f4(x,u1,u2,w1,w2)
output = -u2-u1+*sin(x);

龙格库塔迭代函数

function [u1,u2,u3,w1,w2,w3] = RK4_2variable(u1 , u2 , w1 , w2 , h , a , b , P1 , P2 , T)x = a:h:b;for i = 1:length(x)-1k11 = f1(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
k21 = f2(x(i) , u1(i) , u2(i) , w1(i) , w2(i) , P1 , P2 , T);
L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i) , P1 , P2 , T);k12 = f1(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
k22 = f2(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2 , P1 , P2 , T);
L12 = f3(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2);
L22 = f4(x(i)+h/2 , u1(i)+h*k11/2 , u2(i)+h*k21/2 , w1(i)+h*L11/2, w2(i)+h*L21/2 , P1 , P2 , T);k13 = f1(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
k23 = f2(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2 , P1 , P2 , T);
L13 = f3(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2);
L23 = f4(x(i)+h/2 , u1(i)+h*k12/2 , u2(i)+h*k22/2 , w1(i)+h*L12/2 , w2(i)+h*L22/2 , P1 , P2 , T);k14 = f1(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
k24 = f2(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23 , P1 , P2 , T);
L14 = f3(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23);
L24 = f4(x(i)+h , u1(i)+h*k13 , u2(i)+h*k23 , w1(i)+h*L13 , w2(i)+h*L23 , P1 , P2 , T);u1(i+1) = u1(i) + h/6 * (k11 + 2*k12 + 2*k13 + k14);
u2(i+1) = u2(i) + h/6 * (k21 + 2*k22 + 2*k23 + k24);
w1(i+1) = w1(i) + h/6 * (L11 + 2*L12 + 2*L13 + L14);
w2(i+1) = w2(i) + h/6 * (L21 + 2*L22 + 2*L23 + L24);u3(i) = k11;
w3(i) = L11;
end
end

主函数

% main func
clear;clc;
u1(1) = 0;
u2(1) = 1;
w1(1) = 1;
w2(1) = 0;
h=0.01;
a = 0;b=20;
t = a:h:b;
[shu2,suoyin2]=unique(t2);
t2_=shu2; PU_=PU(suoyin2);
P1=-cos(t); 
P2=-sin(t)-cos(t); 
[u1,u2,u3,w1,w2,w3] = RK4_2variable(u1,u2,w1,w2,h,a,b,P1,P2,t);

运行主函数即可得到近似结果

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

相关文章

Qt重写webrtc的demo peerconnection

整个demo为: 可以选择多个编码方式: cmake_minimum_required(VERSION 3.5)project(untitled LANGUAGES CXX) set(CMAKE_CXX_STANDARD 20) set(CMAKE_INCLUDE_CURRENT_DIR ON)set(CMAKE_AUTOUIC ON) set(CMAKE_AUTOMOC ON) set(CMAKE_AUTORCC ON)set(CMA…

leetcode - 916. Word Subsets

Description You are given two string arrays words1 and words2. A string b is a subset of string a if every letter in b occurs in a including multiplicity. For example, “wrr” is a subset of “warrior” but is not a subset of “world”. A string a from …

如何将 sqlserver 数据迁移到 mysql

文章目录 前言一、导出SQL Server 数据二、转换数据格式为MySQL兼容格式三、导入数据到MySQL数据库五、使用ETL工具六、通过 navicat 工具七、总结 前言 将 SQL Server 数据迁移到 MySQL 是一个常见的数据库迁移任务,通常涉及以下几个关键步骤:导出 SQL…

Pycharm 使用教程

一、基本配置 1. 切换Python解释器 pycharm切换解释器版本 2. pycharm虚拟环境配置 虚拟环境的目的:创建适用于该项目的环境,与系统环境隔离,防止污染系统环境(包括需要的库)虚拟环境配置存放在项目根目录下的 ven…

鸿蒙面试 2025-01-11

ArkTs 和TS的关系? ArkTS(方舟开发语言)与 TypeScript(TS)存在紧密联系,同时也有显著区别: 联系 语法基础:ArkTS 在语法层面大量借鉴了 TypeScript ,TypeScript 里诸如…

蓝桥杯历届真题 #食堂(C++,Java)

这题没什么好说的 考虑所有情况然后写就完了 虽然赛场上 交完不知道答案(doge) 原题链接 #include<iostream>using namespace std;int main() {int n;cin >> n;//能优先安排6人桌,要先安排6人桌//6人桌可以是222 或者 33 或者42//优先用33组合,因为3人寝只能凑6人…

Bash语言的语法糖

Bash语言的语法糖 引言 在现代编程语言中&#xff0c;“语法糖”是一个非常常见的术语&#xff0c;它指的是那些使代码更加易读、易写的语法特性。尽管这些特性并不改变语言的功能&#xff0c;但它们能显著提升开发者的编程体验。在众多编程语言中&#xff0c;Bash&#xff0…

Linux -- 自定义协议体会序列化和反序列化

思路介绍 网络版计算器&#xff1a; 1、客户端发送 两个操作数 和 操作符&#xff1b; 2、根据协议&#xff0c;在发送时&#xff0c;对数据进行序列化&#xff0c;再加上报头&#xff0c;形成 请求 并发送给 服务器&#xff1b; 3、服务器收到 请求 后&#xff0c;判断收到的 …