matlab实践(一):利用ode45和四阶龙哥库塔解二阶耦合微分方程

2023-12-31 21:50

本文主要是介绍matlab实践(一):利用ode45和四阶龙哥库塔解二阶耦合微分方程,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

1.题目

\begin{aligned} \frac{d^2R_1}{dt^2}+\frac32\frac1{R_1}\bigg(\frac{dR_1}{dt}\bigg)^2+\frac{4\mu}{\rho R_1^2}\frac{dR_1}{dt}+\frac{R_2^2}{R_1L}\bigg[\frac{d^2R_2}{dt^2}+\frac2{R_2}\bigg(\frac{dR_2}{dt}\bigg)^2\bigg] \\ =\frac{1}{R_{1}\rho}\Bigg[p_{_{v}}-p_{_{0}}-\frac{2\sigma}{R_{_{1}}}+p_{_{1g0}}\Bigg(\frac{R_{_{10}}}{R_{_{1}}}\Bigg)^{3n}-P\Bigg] \\ \frac{d^{2}R_{2}}{dt^{2}}+\frac{3}{2}\frac{1}{R_{2}}\bigg(\frac{dR_{2}}{dt}\bigg)^{2}+\frac{4\mu}{\rho R_{2}^{2}}\frac{dR_{2}}{dt}+\frac{R_{1}^{2}}{R_{z}L}\bigg[\frac{d^{2}R_{1}}{dt^{z}}+\frac{2}{R_{_1}}\bigg(\frac{dR_{_1}}{dt}\bigg)^{2}\bigg]\\ =\frac{1}{R_{2}\rho}\Bigg[p_{_{v}}-p_{_{0}}-\frac{2\sigma}{R_{_{2}}}+p_{_{2g0}}\Bigg(\frac{R_{_{20}}}{R_{_{2}}}\Bigg)^{3n}-P\Bigg] \\ \end{aligned}

2.ode45

2.1工具箱介绍

ode45 - 求解非刚性微分方程 - 中阶方法

    此 MATLAB 函数(其中 tspan = [t0 tf])求微分方程组 y'=f(t,y) 从 t0 到 tf 的积分,初始条件为 y0。解数组 y
    中的每一行都与列向量 t 中返回的值相对应。

    [t,y] = ode45(odefun,tspan,y0)
    [t,y] = ode45(odefun,tspan,y0,options)
    [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options)
    sol = ode45(___)

2.2求解过程

利用{y_{1}}'=f_{1},{y_{2}}'=f_{2},将这二阶方程化为四个方程,编写函数利用ode45求解,得出结果。

function dy=kongqi(t,y)
n=1.4;
p=1e3;
p_o=2e5;
p_v=1e5;
u=1e-3;
b=7.061e-2;
p_1go=2e4;
p_2go=2e4;
a1=1e-4;
a2=1e-4;
P=5e5*sin(2e4*pi*t);
L=1e-3;
dy=zeros(4,1);
dy(1)=y(2);
dy(2)=-1.5/y(1)*y(2)^2-4*u/(p*y(1)^2)*y(2)-y(3)^2/(y(1)*L)*(dy(4)+2/y(3)*y(4)^2)+1/(y(1)*p)*(p_v-p_o-2*b/y(1)+p_1go*(a1/y(1))^(3*n)-P);
dy(3)=y(4);
dy(4)=-1.5/y(3)*y(4)^2-4*u/(p*y(3)^2)*y(4)-y(1)^2/(y(3)*L)*(dy(2)+2/y(1)*y(2)^2)+1/(y(3)*p)*(p_v-p_o-2*b/y(3)+p_2go*(a2/y(3))^(3*n)-P);
end

这是调用ode45求解。

clc;clear;
tspan=[0 1.5e-5];
y0=[1e-4,0,1e-4,0];
[t1,y1] = ode45('kongqi',tspan,y0);plot(t1,1e4*y1(:,1));

2.3结果

3.四阶龙哥库塔

我们可以利用自己编写的四阶龙哥库塔来求解。

3.1理论知识

经典的四阶龙哥库塔公式如下:

\begin{cases}y_{n+1}=y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4)\\\\K_1=F(x_n,y_n)\\\\K_2=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1)\\\\K_3=f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2)\\\\K_4=F(x_n+h,y_n+hK_3)\end{cases}

其中K1,K2,K3,K4为不同函数值。

3.2代码

这是四阶龙哥库塔函数的代码

function [u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b)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));
L11 = f3(x(i) , u1(i) , u2(i) , w1(i) , w2(i));
L21 = f4(x(i) , u1(i) , u2(i) , w1(i) , w2(i));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);
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);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);
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);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);
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);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);end
end

我们编写四个微分方程:

function output = f1(x,u1,u2,w1,w2)
output = u2;
end
function output = f2(x,u1,u2,w1,w2)
n=1.4;
l=0.1;
P_a=5e5;
R0=1e-4;
A=4e-7*l/R0;
B=7.061e-7*l/R0;
c=1e-5*P_a;
k=0.2;
output =(1.5*l^2*(l*w1*w2^2-u2^2)+2*l^3*(l*u2^2*w1-w2^2)-k*(l*w1^(1-3*n)-u1^(-3*n))+A*(l*w2-u2/u1))/(l^2*u1-l^4*u1^2*w1)+(2*B*(1-1/(l*u1))+(1+c)*(l*w1-1))/(l^2*u1-l^4*u1^2*w1);
end
function output = f3(x,u1,u2,w1,w2)
output = w2;
end
function output=f4(x,u1,u2,w1,w2)
n=1.4;
l=0.1;
P_a=5e5;
R0=1e-4;A=4e-7*l/R0;
B=7.061e-7*l/R0;
c=1e-5*P_a;
k=0.2;
output=(1.5*l^2*(l*u1*u2^2-w2^2)+2*l^3*(l*w2^2*u1-u2^2)-k*(l*u1^(1-3*n)-w1^(-3*n))+A*(l*u2-w2/w1))/(l^2*w1-l^4*w1^2*u1)+(2*B*(1-1/(l*w1))+(1+c)*(l*u1-1))/(l^2*w1-l^4*w1^2*u1);
end

最后利用自己所写的函数求解:

clc;clear;
u1(1) = 1;
u2(1) = 0;
w1(1) = 1;
w2(1) = 0;
h=0.0001;
a = 0;b=0.15;
[u1,u2,w1,w2] = RK4_2variable(u1,u2,w1,w2,h,a,b);
plot(a:h:b,u1,'b-');

3.3结果

这篇关于matlab实践(一):利用ode45和四阶龙哥库塔解二阶耦合微分方程的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



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

相关文章

Spring Boot集成/输出/日志级别控制/持久化开发实践

《SpringBoot集成/输出/日志级别控制/持久化开发实践》SpringBoot默认集成Logback,支持灵活日志级别配置(INFO/DEBUG等),输出包含时间戳、级别、类名等信息,并可通过... 目录一、日志概述1.1、Spring Boot日志简介1.2、日志框架与默认配置1.3、日志的核心作用

破茧 JDBC:MyBatis 在 Spring Boot 中的轻量实践指南

《破茧JDBC:MyBatis在SpringBoot中的轻量实践指南》MyBatis是持久层框架,简化JDBC开发,通过接口+XML/注解实现数据访问,动态代理生成实现类,支持增删改查及参数... 目录一、什么是 MyBATis二、 MyBatis 入门2.1、创建项目2.2、配置数据库连接字符串2.3、入

Android Paging 分页加载库使用实践

《AndroidPaging分页加载库使用实践》AndroidPaging库是Jetpack组件的一部分,它提供了一套完整的解决方案来处理大型数据集的分页加载,本文将深入探讨Paging库... 目录前言一、Paging 库概述二、Paging 3 核心组件1. PagingSource2. Pager3.

在Java中使用OpenCV实践

《在Java中使用OpenCV实践》用户分享了在Java项目中集成OpenCV4.10.0的实践经验,涵盖库简介、Windows安装、依赖配置及灰度图测试,强调其在图像处理领域的多功能性,并计划后续探... 目录前言一 、OpenCV1.简介2.下载与安装3.目录说明二、在Java项目中使用三 、测试1.测

MyBatis-Plus 自动赋值实体字段最佳实践指南

《MyBatis-Plus自动赋值实体字段最佳实践指南》MyBatis-Plus通过@TableField注解与填充策略,实现时间戳、用户信息、逻辑删除等字段的自动填充,减少手动赋值,提升开发效率与... 目录1. MyBATis-Plus 自动赋值概述1.1 适用场景1.2 自动填充的原理1.3 填充策略

Olingo分析和实践之EDM 辅助序列化器详解(最佳实践)

《Olingo分析和实践之EDM辅助序列化器详解(最佳实践)》EDM辅助序列化器是ApacheOlingoOData框架中无需完整EDM模型的智能序列化工具,通过运行时类型推断实现灵活数据转换,适用... 目录概念与定义什么是 EDM 辅助序列化器?核心概念设计目标核心特点1. EDM 信息可选2. 智能类

Olingo分析和实践之OData框架核心组件初始化(关键步骤)

《Olingo分析和实践之OData框架核心组件初始化(关键步骤)》ODataSpringBootService通过初始化OData实例和服务元数据,构建框架核心能力与数据模型结构,实现序列化、URI... 目录概述第一步:OData实例创建1.1 OData.newInstance() 详细分析1.1.1

Olingo分析和实践之ODataImpl详细分析(重要方法详解)

《Olingo分析和实践之ODataImpl详细分析(重要方法详解)》ODataImpl.java是ApacheOlingoOData框架的核心工厂类,负责创建序列化器、反序列化器和处理器等组件,... 目录概述主要职责类结构与继承关系核心功能分析1. 序列化器管理2. 反序列化器管理3. 处理器管理重要方

虚拟机Centos7安装MySQL数据库实践

《虚拟机Centos7安装MySQL数据库实践》用户分享在虚拟机安装MySQL的全过程及常见问题解决方案,包括处理GPG密钥、修改密码策略、配置远程访问权限及防火墙设置,最终通过关闭防火墙和停止Net... 目录安装mysql数据库下载wget命令下载MySQL安装包安装MySQL安装MySQL服务安装完成

SpringBoot整合(ES)ElasticSearch7.8实践

《SpringBoot整合(ES)ElasticSearch7.8实践》本文详细介绍了SpringBoot整合ElasticSearch7.8的教程,涵盖依赖添加、客户端初始化、索引创建与获取、批量插... 目录SpringBoot整合ElasticSearch7.8添加依赖初始化创建SpringBoot项