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的配置文件是灵活且强大的工具,通过合理的配置管理,可以让应用开发和部署更加高效,无论是简单的属性配置,还是复杂... 目录Spring Boot 配置文件详解一、Spring Boot 配置文件类型1.1 applicatio

tomcat多实例部署的项目实践

《tomcat多实例部署的项目实践》Tomcat多实例是指在一台设备上运行多个Tomcat服务,这些Tomcat相互独立,本文主要介绍了tomcat多实例部署的项目实践,具有一定的参考价值,感兴趣的可... 目录1.创建项目目录,测试文China编程件2js.创建实例的安装目录3.准备实例的配置文件4.编辑实例的

Python 中的异步与同步深度解析(实践记录)

《Python中的异步与同步深度解析(实践记录)》在Python编程世界里,异步和同步的概念是理解程序执行流程和性能优化的关键,这篇文章将带你深入了解它们的差异,以及阻塞和非阻塞的特性,同时通过实际... 目录python中的异步与同步:深度解析与实践异步与同步的定义异步同步阻塞与非阻塞的概念阻塞非阻塞同步

Python Dash框架在数据可视化仪表板中的应用与实践记录

《PythonDash框架在数据可视化仪表板中的应用与实践记录》Python的PlotlyDash库提供了一种简便且强大的方式来构建和展示互动式数据仪表板,本篇文章将深入探讨如何使用Dash设计一... 目录python Dash框架在数据可视化仪表板中的应用与实践1. 什么是Plotly Dash?1.1

springboot集成Deepseek4j的项目实践

《springboot集成Deepseek4j的项目实践》本文主要介绍了springboot集成Deepseek4j的项目实践,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价... 目录Deepseek4j快速开始Maven 依js赖基础配置基础使用示例1. 流式返回示例2. 进阶

Android App安装列表获取方法(实践方案)

《AndroidApp安装列表获取方法(实践方案)》文章介绍了Android11及以上版本获取应用列表的方案调整,包括权限配置、白名单配置和action配置三种方式,并提供了相应的Java和Kotl... 目录前言实现方案         方案概述一、 androidManifest 三种配置方式

Spring Boot中定时任务Cron表达式的终极指南最佳实践记录

《SpringBoot中定时任务Cron表达式的终极指南最佳实践记录》本文详细介绍了SpringBoot中定时任务的实现方法,特别是Cron表达式的使用技巧和高级用法,从基础语法到复杂场景,从快速启... 目录一、Cron表达式基础1.1 Cron表达式结构1.2 核心语法规则二、Spring Boot中定

Ubuntu中Nginx虚拟主机设置的项目实践

《Ubuntu中Nginx虚拟主机设置的项目实践》通过配置虚拟主机,可以在同一台服务器上运行多个独立的网站,本文主要介绍了Ubuntu中Nginx虚拟主机设置的项目实践,具有一定的参考价值,感兴趣的可... 目录简介安装 Nginx创建虚拟主机1. 创建网站目录2. 创建默认索引文件3. 配置 Nginx4

Nginx实现高并发的项目实践

《Nginx实现高并发的项目实践》本文主要介绍了Nginx实现高并发的项目实践,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习价值,需要的朋友们下面随着小编来一起学习学习吧... 目录使用最新稳定版本的Nginx合理配置工作进程(workers)配置工作进程连接数(worker_co

Spring Retry 实现乐观锁重试实践记录

《SpringRetry实现乐观锁重试实践记录》本文介绍了在秒杀商品SKU表中使用乐观锁和MybatisPlus配置乐观锁的方法,并分析了测试环境和生产环境的隔离级别对乐观锁的影响,通过简单验证,... 目录一、场景分析 二、简单验证 2.1、可重复读 2.2、读已提交 三、最佳实践 3.1、配置重试模板