如何利用MATLAB代码求解微分方程组

最近在编写代码时遇到了问题，如何求下式的微分方程组。目的是求出振幅-频率的响应曲线，式中w=2πf，f为频率，范围在35-50，t为时间。

其大致参数如下

kn为237.6，a为3*10^(-5)m。
fx为分段函数，代码如下

function fx=fx(t)
d1=2;d2=2.38;d3=2.63;d4=2.83;
e1=0.5;e2=1;e3=1.5;e4=2;
if t>=-e1 && t<e1    
fx=0;
elseif(t>=e1 && t<e2)    
fx=d1^4*(t-e1);
elseif(t>=e2 && t<=e3)    
fx=(d1^4*(t-e1)+d2^4*(t-e2));
elseif(t>e3 && t<=e4)    
fx=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3));
elseif(t>e4)    
fx=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3)+d4^3*(t-e4));
elseif(t>=-e2 && t<-e1)    
fx=-d1^4*(-t-e1);elseif(t>=-e3 && t<-e2)   
 fx=-(d1^4*(-t-e1)+d2^4*(-t-e2));
elseif(t>=-e4 && t<-e3)    fx=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3));
elseif(t<-e4)    
fx=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3)+d4^3*(-t-e4));
end

该回答引用ChatGPT
MATLAB可以通过ode函数求解微分方程组，其语法为：

[T,Y] = ode45(ODEFUN,TSPAN,Y0)

其中，

• T 是时间向量
• Y 是解向量，它的每一行都给出一个与相应时间处的解。
• ODEFUN 指的是调用函数，它需要接受两个输入变量：t 和 y。t 是当前时间，y 是对应的状态向量。
• TSPAN 是一个包含求解时间的向量，例如：TSPAN=linspace(t0,t1,tsteps)，其中 t0 和 t1 是开始和结束时间，tsteps 是时间步数。
• Y0 是初始条件向量。

下面是一个求解微分方程组的简单例子：

function dydt = myODE(t,y)
dydt = zeros(2,1);
dydt(1) = -20*y(1) + 40*y(2) + 30*sin(t);
dydt(2) = 20*y(1) - 40*y(2) + 10*cos(t);

[t,y] = ode45(@myODE,[0 10],[1 0]);

在这个例子中，myODE 函数返回一个列向量 dydt，其中 dydt(1) 和 dydt(2) 分别对应第一个和第二个状态变量的导数。ode45 函数用于求解微分方程组，它传递给 myODE 函数的时间向量为 t，与状态向量为 y。最后，例如对于时间从0到10的时间步的初值问题，y的初值为 [1 0]。

参考GPT：要求解微分方程组，可以使用MATLAB中的ode45函数。这个函数可以求解常微分方程组（ODE），并返回一个包含解的矩阵的结构。下面是一个示例代码，用于求解您提供的微分方程组：

% 定义参数
m = 1;
mn = 0.08;
cx = 6;
cy = 6;
kx = 76400;
ky = 76400;
ko = 2000;
kn = 237.6;
a = 3e-5;
 
% 定义初始条件
x0 = 0;
v0x = 0;
y0 = 0;
v0y = 0;
xn0 = 0;
v0xn = 0;
yn0 = 0;
v0yn = 0;
 
% 定义ODE函数
odefun = @(t, y) [
    y(2); % dx/dt
    (m*a^2*cos(2*pi*t) - cx*y(2) - kx*y(1) - ko*(y(1) - y(5)))/m; % d2x/dt2
    y(4); % dy/dt
    (m*a^2*sin(2*pi*t) - cy*y(4) - ky*y(3) - ko*(y(3) - y(6)))/m; % d2y/dt2
    y(6); % dxn/dt
    (-kn*y(5) - ko*(y(5) - y(1)) - cx*y(6) + fx(t))/mn; % d2xn/dt2
    y(8); % dyn/dt
    (-kn*y(7) - ko*(y(7) - y(3)) - cy*y(8) + fy(t))/mn; % d2yn/dt2
];
 
% 求解ODE
tspan = [0 1]; % 求解时间范围
[t, sol] = ode45(odefun, tspan, [x0 v0x y0 v0y xn0 v0xn yn0 v0yn]);
 
% 绘制响应曲线
x = sol(:, 1);
y = sol(:, 3);
plot(x, y);
xlabel('振幅');
ylabel('频率');
function fx=fx(t)
d1=2;d2=2.38;d3=2.63;d4=2.83;
e1=0.5;e2=1;e3=1.5;e4=2;
if t>=-e1 && t<e1    
fx=0;
elseif(t>=e1 && t<e2)    
fx=d1^4*(t-e1);
elseif(t>=e2 && t<=e3)    
fx=(d1^4*(t-e1)+d2^4*(t-e2));
elseif(t>e3 && t<=e4)    
fx=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3));
elseif(t>e4)    
fx=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3)+d4^3*(t-e4));
elseif(t>=-e2 && t<-e1)    
fx=-d1^4*(-t-e1);elseif(t>=-e3 && t<-e2)   
 fx=-(d1^4*(-t-e1)+d2^4*(-t-e2));
elseif(t>=-e4 && t<-e3)    fx=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3));
elseif(t<-e4)    
fx=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3)+d4^3*(-t-e4));
end
end
function fy=fy(t)
d1=2;d2=2.38;d3=2.63;d4=2.83;
e1=0.5;e2=1;e3=1.5;e4=2;
if t>=-e1 && t<e1    
fy=0;
elseif(t>=e1 && t<e2)    
fy=d1^4*(t-e1);
elseif(t>=e2 && t<=e3)    
fy=(d1^4*(t-e1)+d2^4*(t-e2));
elseif(t>e3 && t<=e4)    
fy=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3));
elseif(t>e4)    
fy=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3)+d4^3*(t-e4));
elseif(t>=-e2 && t<-e1)    
fy=-d1^4*(-t-e1);elseif(t>=-e3 && t<-e2)   
 fy=-(d1^4*(-t-e1)+d2^4*(-t-e2));
elseif(t>=-e4 && t<-e3)    fy=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3));
elseif(t<-e4)    
fy=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3)+d4^3*(-t-e4));
end
end

其中，ODE函数中的y向量是微分方程组的解，第一个元素是x，第二个元素是dx/dt，第三个元素是y，第四个元素是dy/dt，第五个元素是xn，第六个元素是dxn/dt，第七个元素是yn，第八个元素是dyn/dt。fx(t)和fy(t)是您提供的分段函数fx和fy在时间t处的值。

在求解ODE后，您可以从解矩阵sol中提取x和y向量，并绘制它们的响应曲线。

该回答引用ChatGPT GPT-4

运行结果

函数

function dydt = sys_diffeq(t, y, m, ma, cx, cy, k, ko, fx, fy, kn, mn, cn, w)

    dydt = zeros(8, 1);

    dydt(1) = y(2);
    dydt(2) = (ma*w^2*cos(w*t) - cx*y(2) - k*y(1) - ko*(y(1)-y(5)))/m;
    dydt(3) = y(4);
    dydt(4) = (ma*w^2*sin(w*t) - cy*y(4) - k*y(3) - ko*(y(3)-y(7)))/m;
    dydt(5) = y(6);
    dydt(6) = (-cn*y(6) - ko*(y(5)-y(1)) - fx*(y(5)-y(1)) - kn*y(5))/mn;
    dydt(7) = y(8);
    dydt(8) = (-cn*y(7) - ko*(y(7)-y(3)) - fy*(y(7)-y(3)) - kn*y(7))/mn;
end

主程序

% 参数
m = 1;
ma = 0.08;
cx = 6;
cy = 3;
k = 76400;
ko = 2000;
fx = 0;
fy = 0;
kn = 0;
mn = 1;
cn = 3;


% 频率范围
f_min = 35;
f_max = 50;
num_freqs = 100;
freqs = linspace(f_min, f_max, num_freqs);

% 时间范围
t_span = [0, 10];

% 初始条件
y0 = [0; 0; 0; 0; 0; 0; 0; 0];

% 振幅-频率曲线
amplitudes = zeros(num_freqs, 1);

for i = 1:num_freqs
    f = freqs(i);
    w = 2 * pi * f;

    % 使用ode45求解微分方程组
   [t, y] = ode45(@(t, y) sys_diffeq(t, y, m, ma, cx, cy, k, ko, fx, fy, kn, mn, cn, w), t_span, y0);


    % 计算振幅
    x_amplitude = max(y(:, 1)) - min(y(:, 1));
    y_amplitude = max(y(:, 3)) - min(y(:, 3));
    amplitudes(i) = sqrt(x_amplitude^2 + y_amplitude^2);
end

figure;
plot(freqs, amplitudes);
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Amplitude-Frequency Response');
grid on;

• 可以看下matlab参考手册中的 matlab 求解刚性微分方程 - 低阶方法 ode23s

该回答引用于gpt与OKX安生共同编写：

• 该回答引用于gpt与OKX安生共同编写：

您可以使用 MATLAB 中的 ode45 函数来求解微分方程组。在使用 ode45 前，需要将微分方程组转化为向量形式。

对于您给出的微分方程组，可以将其转化为以下形式：

function dydt = myode(t,y)
    % y(1) = x
    % y(2) = v
    kn = 237.6;
    a = 3e-5;
    fx = fx(t);
    w = 2 * pi * linspace(35,50,100); % 生成频率范围数组
    F = fft(fx, 100) / 100; % 对力函数进行 FFT
    F = F(1:length(F)/2); % 取前一半的幅度谱
    A = zeros(size(w)); % 频率响应的振幅
    for i = 1:length(w)
        H = (kn/m) ./ (-w(i)^2 + 1i*w(i)*b + k/m); % 计算传递函数
        A(i) = abs(H*F(i)); % 计算振幅
    end
    dydt = [y(2); -k/m*y(1) - b/m*y(2) + fx(t)/m]; % 状态变量的微分
end

其中，dydt 是状态变量 x 和 v 的微分向量。h 是传递函数，A 是频率响应的振幅。

然后，在主程序中调用 ode45 函数求解微分方程组：

tspan = [0, 10]; % 时间范围
y0 = [0, 0]; % 初始状态
[t, y] = ode45(@myode, tspan, y0); % 求解微分方程组
plot(t, y(:,1)); % 画出位移随时间变化的图像
xlabel('Time (s)');
ylabel('Displacement (m)');

在这里，我们将初始状态设为 $x(0)=0$ 和 $v(0)=0$。然后，我们使用 ode45 函数求解微分方程组，并用 plot 函数画出模拟结果。

• 如这个实现对您有所帮助的话，还请采纳哦~

要求求解微分方程组，可以使用MATLAB中的ode45函数进行求解。首先需要将微分方程组转化为MATLAB可以处理的形式。
将式子中的x1、x2、x3、x4分别表示为y(1)、y(2)、y(3)、y(4)，则微分方程组可以表示为：
dy(1)/dt = y(2)
dy(2)/dt = -2*zeta*wn*y(2) - wn^2*y(1) + fx(t)/m
dy(3)/dt = y(4)
dy(4)/dt = -2*zeta*wn*y(4) - wn^2*y(3)
其中，wn为自然频率，zeta为阻尼比，m为质量。
然后，可以编写MATLAB代码进行求解。代码如下：

matlab
function dydt = myode(t,y)
% 参数设置
wn = 2*pi*linspace(35,50,100); % 自然频率
zeta = 0.01; % 阻尼比
m = 1; % 质量
% 分段函数
function fx = fx(t)
d1=2;d2=2.38;d3=2.63;d4=2.83;
e1=0.5;e2=1;e3=1.5;e4=2;
if t>=-e1 && t    fx=0;
elseif(t>=e1 && t    fx=d1^4*(t-e1);
elseif(t>=e2 && t<=e3)
    fx=(d1^4*(t-e1)+d2^4*(t-e2));
elseif(t>e3 && t<=e4)
    fx=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3));
elseif(t>e4)
    fx=(d1^4*(t-e1)+d2^4*(t-e2)+d3^4*(t-e3)+d4^3*(t-e4));
elseif(t>=-e2 && t<-e1)
    fx=-d1^4*(-t-e1);
elseif(t>=-e3 && t<-e2)
    fx=-(d1^4*(-t-e1)+d2^4*(-t-e2));
elseif(t>=-e4 && t<-e3)
    fx=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3));
elseif(t<-e4)
    fx=-(d1^4*(-t-e1)+d2^4*(-t-e2)+d3^4*(-t-e3)+d4^3*(-t-e4));
end
% 求解微分方程组
dydt = zeros(4,1);
dydt(1) = y(2);
dydt(2) = -2*zeta*wn*y(2) - wn.^2*y(1) + fx(t)/m;
dydt(3) = y(4);
dydt(4) = -2*zeta*wn*y(4) - wn.^2*y(3);
end
% 初始条件
y0 = [0;0;0;0];
% 求解微分方程组
[t,y] = ode45(@myode,[0,10],y0);
% 绘制振幅-频率响应曲线
A = sqrt(y(:,1).^2 + y(:,3).^2);
plot(wn/(2*pi),A);
xlabel('频率（Hz）');
ylabel('振幅（m）');