matlab实践（一）：利用ode45和四阶龙哥库塔解二阶耦合微分方程

您所在的位置：网站首页 › matlab解二阶微分 › matlab实践（一）：利用ode45和四阶龙哥库塔解二阶耦合微分方程

matlab实践（一）：利用ode45和四阶龙哥库塔解二阶耦合微分方程

2024-07-09 10:11:44| 来源: 网络整理| 查看: 265

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)-1 k11 = 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结果

【本文地址】

公司简介

联系我们

今日新闻

点击排行

实验室常用的仪器、试剂和: 说到实验室常用到的东西，主要就分为仪器、试剂和耗

不用再找了，全球10大实验: 01、赛默飞世尔科技（热电）Thermo Fisher Scientif

三代水柜的量产巅峰T-72坦: 作者：寞寒最近，西边闹腾挺大，本来小寞以为忙完这

通风柜跟实验室通风系统有: 说到通风柜跟实验室通风，不少人都纠结二者到底是不

集消毒杀菌、烘干收纳为一: 厨房是家里细菌较多的地方，潮湿的环境、没有完全密

实验室设备之全钢实验台如: 全钢实验台是实验室家具中较为重要的家具之一，很多

图片新闻

实验室药品柜的特性有哪些: 实验室药品柜是实验室家具的重要组成部分之一，主要

小学科学实验中有哪些教学: 计算机计算器一般打孔器打气筒仪器车显微镜

实验室各种仪器原理动图讲: 1.紫外分光光谱UV分析原理：吸收紫外光能量，引起分

高中化学常见仪器及实验装: 1、可加热仪器：2、计量仪器：（1）仪器A的名称：量

微生物操作主要设备和器具: 今天盘点一下微生物操作主要设备和器具，别嫌我啰嗦

浅谈通风柜使用基本常识: 　众所周知，通风柜功能中最主要的就是排气功能。在

matlab实践（一）：利用ode45和四阶龙哥库塔解二阶耦合微分方程

matlab实践（一）：利用ode45和四阶龙哥库塔解二阶耦合微分方程

今日新闻

点击排行

推荐新闻

图片新闻

专题文章