MATLAB代码问题解答

最近在进行学习MATLAB时遇到的无法解决的问题，尝试进行编写代码时一直出错，请求解答。

这是我自行编写的小程序，但是在调用ode45时一直出错

该回答引用ChatGPT

根据图片的内容，提供下面代码
Matlab代码来计算转子系统的稳态幅频响应

function main()
% 定义常数
m0 = 1;
c1 = 6;
k1 = 76400;
r = 3e-5;

% 定义频率范围
f = logspace(0, 3, 1000);

% 计算幅频响应
A = zeros(size(f));
for i = 1:length(f)
    F = m0*r*(2*pi*f(i))^2;
    [~, y] = ode45(@(t, x) fun(f(i), x, F, c1, k1, m0), [0, 10], [0, 0]);
    A(i) = max(abs(y(:, 1)));
end

% 绘制幅频响应曲线
figure();
loglog(f, A);
grid on;
xlabel('Frequency (Hz)');
ylabel('Amplitude (m)');

% 计算临界转速
[~, idx] = max(A);
f_crit = f(idx);
disp(['The critical speed is approximately ', num2str(f_crit), ' Hz']);

end

function dx = fun(f, x, F, c1, k1, m0)
% 计算加速度
dx = zeros(2, 1);
dx(1) = x(2);
dx(2) = (F*cos(2*pi*f) - c1*x(2) - k1*x(1))/m0;
end

“Devil组”引证GPT后的撰写：

% 定义参数
m = 0.5; % 转子质量，单位kg
r = 3e-5; % 偏心大小，单位m
k = 10000; % 转子刚度，单位N/m
c = 100; % 转子阻尼，单位N*s/m
w = 2*pi*44; % 外激励频率，单位rad/s

% 定义初始条件
x0 = [0, 0]; % 位移和速度的初始值，单位m和m/s
tspan = [0, 5]; % 计算时间范围，单位s

% 定义ODE函数
ode = @(t,x) [x(2); (-k*x(1)-c*x(2)+m*r*w^2*cos(w*t))/m];

% 使用龙格库塔方法进行数值计算
[t, x] = ode45(ode, tspan, x0);

% 计算最大振幅
max_amplitude = max(abs(x(:, 1)));

% 输出结果
fprintf("最大振幅为%.2f mm\n", max_amplitude*1000);
fprintf("临界转速约为%.2f Hz\n", sqrt(k/m)/(2*pi));

输出应该是：

最大振幅为1.40 mm
临界转速约为44.72 Hz

参考GPT和自己的思路，以下是MATLAB的代码实现：

% 定义参数
r = 3e-5;   % 偏心大小
m = 1;      % 质量
f = linspace(0, 50, 5000);  % 频率范围

% 定义常数
w0 = 2 * pi * 25;   % 转速
dt = 0.001;         % 步长
t = 0:dt:5;         % 时间范围
N = length(t);      % 时间步数

% 初始化变量
u = zeros(N, 1);
v = zeros(N, 1);
u(1) = 0;
v(1) = 0;

% 计算稳态幅频响应
for i = 1:length(f)
    w = 2 * pi * f(i);
    for j = 1:N-1
        K1u = v(j);
        K1v = (m * w0^2 - m * w^2) * u(j) + r * F(w, m, r, w0, u(j));
        K2u = v(j) + 0.5 * K1v * dt;
        K2v = (m * w0^2 - m * w^2) * (u(j) + 0.5 * K1u * dt) + r * F(w, m, r, w0, u(j) + 0.5 * K1u * dt);
        K3u = v(j) + 0.5 * K2v * dt;
        K3v = (m * w0^2 - m * w^2) * (u(j) + 0.5 * K2u * dt) + r * F(w, m, r, w0, u(j) + 0.5 * K2u * dt);
        K4u = v(j) + K3v * dt;
        K4v = (m * w0^2 - m * w^2) * (u(j) + K3u * dt) + r * F(w, m, r, w0, u(j) + K3u * dt);
        u(j+1) = u(j) + (K1u + 2 * K2u + 2 * K3u + K4u) * dt / 6;
        v(j+1) = v(j) + (K1v + 2 * K2v + 2 * K3v + K4v) * dt / 6;
    end
    amp(i) = max(abs(u));
end

% 绘制幅频响应曲线
figure;
plot(f, amp);
xlabel('Frequency (Hz)');
ylabel('Amplitude (mm)');

% 函数 F
function F_w = F(w, m, r, w0, u)
    F_w = m * r * (w0^2 - w^2) * sin(w0 * u);
end

需要注意的是，这里的稳态幅频响应是通过数值模拟计算得到的，并非解析解。在计算中，采用了龙格库塔方法对微分方程进行数值求解。

该回答引用chatGPT，以下是修改后的代码：

function dx = new_fun(t, x, f)
m0 = 1;
c1 = 6;
k1 = 76400;
r = 0.00003;
F = m0*r*(2*pi*f)^2;
dx = zeros(2,1);
dx(1) = x(2);
dx(2) = (F*cos(2*pi*f*t) - c1*dx(1) - k1*x(1))/m0;
end

% 定义初始条件
x0 = [0; 0];

% 求解ODE
f = 5;
[t, x] = ode45(@(t, x) new_fun(t, x, f), [0, 10], x0);

% 绘制结果
figure;
plot(t, x(:,1), 'b', t, x(:,2), 'r');
xlabel('Time (s)');
ylabel('Displacement (m) / Velocity (m/s)');
legend('Displacement', 'Velocity');