利用Matlab求解高阶微分方程（ode45）

1、高阶微分方程的基本概念

二阶以及二阶以上的微分方程称之为高阶微分方程，一般来说，微分方程的阶数越高，求解的难度也就越大。求高阶方程的一个常用方法就是降低阶数。对二阶方程，如果能用变量代换把它化成一阶方程，那么就可以用一阶微分方程的求解方法来解决了。

同样，在matlab中不能直接求解高阶微分方程，我们需要使用变量替换将其转换为一阶微分方程来求解。本文基于ode45求解器进行求解。

2、利用ode45求解高阶微分方程

2.1 方法简介

ode45 是 MATLAB 中用于求解常微分方程（ODE）的函数之一，它基于变步长的龙格-库塔方法（5(4)阶方法）。ode45 通过自动调整步长来提高效率和精度。

2.2 函数格式

ode45 的基本语法如下： [t, y] = ode45(odefun, tspan, y0)

odefun: 函数句柄，表示微分方程的右侧函数，应该以 t 和 y 为输入，返回导数 dy/dt。
tspan: 时间范围，通常是一个二元素向量 [t0, tf]，表示时间区间从 t0 到 tf。
y0: 初始条件，是一个向量，包含了在 t0 时刻的解的初始值。
t: 求解出的时间向量，表示在每个时间点上计算的结果。
y: 对应的解向量，与 t 的大小相同。

2.3 举例分析

以2022年国赛A题波浪能最大输出功率设计第一问为例（有些公式出现的很突兀，建议看一下论文，更易理解）参考资料：2022高教社杯全国大学生数学建模竞赛A题论文展示（A001）

第一问建立的数学模型如下：可以发现这是一个二阶二元微分方程， $x_{f},x_{z}$ 是未知量

我们对上述方程进行变化，化为一阶四元微分方程

下面开始编程实现：

编写函数句柄：odefun

function dydt = odefun(t, y, params)% 提取参数zf = params.zf;w = params.w;K_zn = params.K_zn;K_tang = params.K_tang;C_x = params.C_x;B = params.B;S = params.S;mzz = params.mzz;% 提取状态变量x_f = y(1); % 浮子位移x_z = y(2); % 振子位移u = y(3);   % 浮子速度w_n = y(4); % 振子速度% 定义微分方程dx_f = u;dx_z = w_n;du = (zf*cos(w*t) + K_zn*(w_n - u) + K_tang*(x_z - x_f) - C_x*u - B*x_f) / S;dw = -(K_zn*(w_n - u) + K_tang*(x_z - x_f)) / mzz;% 返回结果dydt = [dx_f; dx_z; du; dw];
end

2. 调用ode45函数求解

% 初始化参数
clc; clear; close all;% 系统参数
params.zf = 6250;
params.w = 1.4005;
params.K_zn = 10000;
params.K_tang = 80000;
params.C_x = 656.3616;
params.B = 1025 * 9.8 * pi;
params.S = 4866 + 1335.535;
params.mzz = 2433;% 初始条件 [x_f(0), x_z(0), u(0), w(0)]
y0 = [0, 0, 0, 0];% 时间区间
T = (2 * pi) / params.w * 40; % 前40个波浪周期
tspan = [0 T];% 调用 ode45 求解
[t, y] = ode45(@(t, y) odefun(t, y, params), tspan, y0);% 提取解
x_f = y(:, 1); % 浮子位移
x_z = y(:, 2); % 振子位移
u = y(:, 3);   % 浮子速度
w_n = y(:, 4); % 振子速度% 绘图
figure;% 绘制位移图
subplot(1, 2, 1);
plot(t, x_f, 'r', t, x_z - x_f, 'b');
legend('浮子位移', '振子相对浮子位移');
xlabel('时间 (s)');
ylabel('位移 (m)');
title('位移随时间变化');% 绘制速度图
subplot(1, 2, 2);
plot(t, u, 'r', t, w_n - u, 'b');
legend('浮子速度', '振子相对浮子速度');
xlabel('时间 (s)');
ylabel('速度 (m/s)');
title('速度随时间变化');% 保存结果到 Excel 文件
% 初始化变量
resfx = [];
resfv = [];
reszx = [];
reszv = [];
dt = 0.2; % 时间间隔
time_points = 0:dt:T;% 查找时间点在时间向量中的索引并记录数据
for t_point = time_points[~, idx] = min(abs(t - t_point)); % 查找最接近的时间点resfx(end+1) = x_f(idx);resfv(end+1) = u(idx);reszx(end+1) = x_z(idx);reszv(end+1) = w_n(idx);
end% 保存结果到 Excel 文件
result = [time_points', resfx', resfv', reszx', reszv'];
xlswrite('myresult1.xlsx', result);

结果绘图：