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

1.题目

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求解过程

利用,将这二阶方程化为四个方程,编写函数利用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理论知识

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

其中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结果