超级会员免费看
常见动力学系统的MATLAB建模与分析
1. 霍奇金 - 赫胥黎和菲茨休 - 纳古莫神经元模型
1.1 FHN模型相平面极限环绘制
FHN（FitzHugh - Nagumo）模型是一种用于描述神经元放电行为的数学模型。通过绘制其相平面的极限环，可以直观地观察到模型在不同激励电流下的动态行为。
操作步骤
- 设定激励电流范围：从 0.4 到 1.3，步长为 0.3。
- 定义系统的微分方程。
-
使用
ode45求解器求解微分方程。 - 绘制相平面轨迹。
MATLAB代码
clc;clear;close;syms y(t)
for I=0.4:0.3:1.3
s = @(t,y) [3*(y(1) + y(2) - y(1)^3/3 - I);
(-1/3)*(y(1) - 0.7 + 0.8*y(2))];
[t, Y] = ode45(s, [0 40], [1 0.7]);
plot(Y(:,1),Y(:,2),'linewidth',2.5);hold on;
xlabel('v(t)','FontSize',11);ylabel('w(t)','FontSize',11);grid on;
end
1.2 FHN模型极限环起始与时域曲线绘制
为了研究 FHN 模型极限环的起始条件，需要确定激励电流的上下限 (I_1) 和 (I_2)，并绘制相应的相平面轨迹和时域曲线。
操作步骤
- 设定激励电流范围：从 0.33696 到 0.33698，步长为 0.00001。
- 定义系统的微分方程。
-
使用
ode45求解器求解微分方程。 - 绘制相平面轨迹和时域曲线。
MATLAB代码
clc;clear;close;syms y(t)
for I=0.33696:0.00001:0.33698
s = @(t,y) [3*(y(1) + y(2) - y(1)^3/3 - I);
(-1/3)*(y(1) - 0.7 + 0.8*y(2))];
[t, Y] = ode45(s, [0 40], [0.5 0.2]);
plot(Y(:,1),Y(:,2),'linewidth',2.5);grid; % Limit cycle plot
title(['I = ' num2str(I)]);
xlabel('v(t)','FontSize',11);ylabel('w(t)','FontSize',11);figure;
plot(t,Y(:,1),t,Y(:,2),'k','linewidth',2.5);grid; % v(t) and w(t) vs t
title(['I = ' num2str(I)]);xlabel('t','FontSize',11);
ylabel('v(t), w(t)','FontSize',11);legend('v(t)','w(t)');figure
end
2. 混合槽与间歇反应器中的化学反应
2.1 混合槽中糖质量的动态变化
考虑一个搅拌良好的混合槽，初始时含有一定量的糖和水。在时间为 0 时，以一定浓度和流量的糖溶液连续流入槽中，同时槽内溶液以一定流量流出。
2.1.1 定义系统的常微分方程
通过对槽内糖的质量进行平衡分析，可以得到以下常微分方程：
[
\frac{dM}{dt} = v_{in}C_{in} - v_{out}\frac{M}{V_0 + v_{in}t - v_{out}t}
]
其中，(M) 是槽内糖的质量，(v_{in}) 是流入流量，(C_{in}) 是流入溶液的糖浓度，(v_{out}) 是流出流量，(V_0) 是初始水的体积。
2.1.2 使用
ODE45
求解器绘制糖质量随时间的变化
MATLAB代码
% Numerical solution using ODE45 solver
clc;clear;
tspan = [0,300]; % minutes
V0 = 1000; % L
Mo = 600; % initial sugar mass, g
vin = 8; % inflow (=outflow) rate, L/min
vout = 11; % outflow rate, L/min
Cs = 2; % inflow sugar concentration, g/L
f=@(t,M) vin*Cs - vout*M/(V0+vin*t-vout*t);
[t,M] = ode45(f,tspan,Mo);
plot(t,M,'linewidth',2.5);hold on
xlabel('Time (min)');ylabel('Sugar mass (g)');grid
legend('vout=8L/min','vout=11L/min')
2.1.3 解析解
使用
dsolve
函数进行解析求解：
MATLAB代码
% Analytical solution
clc;clear;
vin=8; Cs=2; V0=1000; syms M(t)
vout=11;
M = dsolve(diff(M)== vin*Cs - vout*M/(V0+vin*t-vout*t), M(0)==600);
M= simplify (M)
h=ezplot(M,[0 , 300]);ylim([0 2000]);
set(h,'color',[0, .45 .74],'linewidth',2.5);hold on;
xlabel('Time (min)');ylabel('Sugar concentration (g/L)');grid on;
legend('vout=8L/min','vout=11L/min');
2.2 间歇反应器中的化学反应动力学优化
在化学生产中，间歇反应器的目标是实现反应过程的最大产品收率。考虑一个假设的双分子反应：
[
a + b \xrightarrow{k_1} c + d
]
[
c + a \xrightarrow{k_2} e
]
操作步骤
- 定义反应速率方程。
- 设定初始浓度。
- 使用数值 ODE 求解器进行模拟。
- 优化反应速率常数 (k_1) 和 (k_2)，以实现产品 (e) 的最大收率。
MATLAB代码
%ode5.m (Chemical) bach reactor
tic;clc;clear;
tf=20;
ic=[1 1 0 0 0];%initial conditions [Ca(0) Cb(0) Cc(0) Cd(0) Ce(0)]
tspan=[0 tf];
k=0;
delta=0.05;
for k1=0.1:delta:1
for k2=0.1:delta:1;
k=k+1;
f = @(t,x) [-k1*x(1)*x(2)-k2*x(1)*x(3);
-k1*x(1)*x(2);
k1*x(1)*x(2)-k2*x(1)*x(3);
k1*x(1)*x(2);
k2*x(1)*x(3)];
[t, Y]=ode23(f,tspan, ic);
a = Y(:,1); b = Y(:,2);c = Y(:,3);d = Y(:,4); e = Y(:,5);
A(k,:)=[k1 k2 Y(end,5)];
end
end
legend('Ca','Cb','Cc','Cd','Ce');grid; %A
[M,I] = max(A(:,3)); disp(' k1, k2, Ce'); A(I,:)
toc
3. 四旋翼动力学建模
3.1 四旋翼的特点与简化模型
四旋翼是一种具有四个垂直方向螺旋桨的旋翼飞行器。与直升机不同，它不需要尾桨，通过改变螺旋桨的推力来实现控制。简化的四旋翼模型可以用一组常微分方程来描述。
3.2 四旋翼模型的常微分方程
四旋翼的动力学模型可以用以下常微分方程描述：
[
\ddot{x} = (\cos\phi \sin\theta \cos\psi + \sin\phi \sin\psi)\frac{1}{m}U_1
]
[
\ddot{y} = (\cos\phi \sin\theta \cos\psi - \sin\phi \cos\psi)\frac{1}{m}U_1
]
[
\ddot{z} = (\cos\phi \cos\theta)\frac{1}{m}U_1 - g
]
[
\ddot{\phi} = \dot{\theta}\dot{\psi}\frac{I_{yy} - I_{zz}}{I_{xx}} + \frac{U_2}{I_{xx}}
]
[
\ddot{\theta} = \dot{\phi}\dot{\psi}\frac{I_{zz} - I_{xx}}{I_{yy}} + \frac{U_3}{I_{yy}}
]
[
\ddot{\psi} = \dot{\phi}\dot{\theta}\frac{I_{xx} - I_{yy}}{I_{zz}} + \frac{U_4}{I_{zz}}
]
其中，(U_1) 是总推力，(U_2) 是俯仰力矩，(U_3) 是滚转力矩，(U_4) 是偏航力矩。
3.3 四旋翼模型的求解与分析
3.3.1 将二阶方程转化为一阶常微分方程组
3.3.2 求解起飞频率
通过令 (z) 方向的加速度为 0，可以得到起飞频率的计算公式：
[
f = \frac{30}{\pi}\sqrt{\frac{mg}{4b}}
]
其中，(m) 是四旋翼的质量，(g) 是重力加速度，(b) 是推力因子。
3.3.3 模拟四旋翼的运动
MATLAB代码
%quadrotor.m Simplified model of a quadrotor (drone)
clc;clear;close;tic
Ixx = 8.1*10^(-3);
Iyy = 8.1*10^(-3);
Izz = 14.2*10^(-3);
b = 54.2*10^(-6);
d = 1.1*10^(-6);
L = 0.24;
m = 1;
g = 9.81;
fc=(30/pi)*(m*g/4/b).^(1/2) % lift-off frequency of each motor
N=40; % number of time segments within given interval
tspan = linspace(0,5,N); % simulation time
x0=zeros(12,1);
f1=2031.4; f2=2031.4; f3=2031.4; f4=2031.4; % revolutions/minute
omega1=2*pi*f1/60; omega2 =2*pi*f2/60;%rad/sec
omega3=2*pi*f3/60; omega4 =2*pi*f4/60;%rad/sec
U1=b*(omega1.^2+omega2.^2+omega3.^2+omega4.^2)
U2=b*L*(-omega1.^2+omega4.^2)
U3=b*L*(omega1.^2-omega3.^2)
U4=d*(-omega1.^2+omega2.^2-omega3.^2+omega4.^2)
% Bounding U
if U1 > 15.7, U1 = 15.7; end;
if U1 < 0, U1 = 0; end;
if U2 > 1, U2 = 1; end;
if U2 < -1, U2 = -1; end;
if U3 > 1 , U3 = 1; end;
if U3 < -1, U3 = -1; end;
if U4 > 1, U4 = 1; end;
if U4 < -1, U4 = -1; end;
fn=@(t,x) [ x(2); % Xdot
(sin(x(11))*sin(x(7)) + cos(x(11))*sin(x(9))*cos(x(7)))*(U1/m);
x(4);
(-cos(x(11))*sin(x(7)) + cos(x(11))*sin(x(9))*cos(x(7)))*(U1/m);
x(6);
(cos(x(9))*cos(x(7)))*(U1/m)-g;
x(8);
((Iyy - Izz)/Ixx)*x(10)*x(12) + (U2/Ixx);
x(10);
((Izz - Ixx)/Iyy)*x(8)*x(12) + (U3/Iyy);
x(12);
((Ixx - Iyy)/Izz)*x(8)*x(10) + (U4/Izz)];
[t,x] = ode45(fn, tspan, x0)
for i=1:N
if x(i,5) < 0, x(i,5) = 0; end
end;
x(:,5)
plot(t,x(:,5),'LineWidth',2.5);% Plot Z-axis motion
title('Altitude'); xlabel('seconds'); ylabel('meters'); grid; figure;
% Plot angle response in degrees
Phid=x(:,7)*180/pi; Thetad=x(:,9)*180/pi; Psid=x(:,11)*180/pi;
plot(t,Phid,'LineWidth',2.5);hold on;
plot(t,Thetad,'r','LineWidth',2.5);hold on;
plot(t,Psid,'k','LineWidth',2.5)
legend('Roll angle ','Pitch angle ','Yaw angle ');
grid; title('Angle response'); xlabel('seconds'); ylabel('Degrees');toc
4. 摆问题
4.1 无阻尼单摆
无阻尼单摆的运动方程为：
[
\frac{d^2\theta}{dt^2} = -\omega^2\sin\theta
]
其中，(\omega = \sqrt{\frac{g}{L}})，(g) 是重力加速度，(L) 是摆长。
操作步骤
-
将二阶方程转化为一阶常微分方程组：
[
\begin{cases}
y_1’ = y_2 \
y_2’ = -\omega^2\sin y_1
\end{cases}
] -
使用
ode45求解器求解微分方程。 - 绘制摆角随时间的变化曲线。
MATLAB代码
%pendulum1.m Solves y"+w^2*sin(y), y=theta, y0=10 degrees, y'(0)=0.
clc;clear;close;
g=9.81;L=1;%SI standard units
w=sqrt(g/L)
f=w/2/pi
T=1/f
tspan=[0 2]; y0=[85 0];%degrees
y0=(pi/180)*y0%radians
pendulum= @(t,y) [y(2); -w^2*sin(y(1))];
[t,y] = ode45(pendulum,tspan, y0);
y=y*180/pi;%degrees
plot(t,y(:,1),'linewidth',2.5);
xlabel('t'); ylabel('y_{1}(t), degrees');title('\theta (t)');grid
figure;
plot(y(:,1),y(:,2),'linewidth',2); xlabel('\theta (t)');grid on;
ylabel('d \theta / dt (t)');
4.2 阻尼驱动摆
阻尼驱动摆的运动方程为：
[
\theta’’ = -\frac{b}{mL^2}\theta’ - \frac{g}{L}\sin\theta + \frac{N_d}{mL^2}\cos\omega_dt
]
通过引入无量纲参数，可以将方程简化为：
[
x’’ = -cx’ - \sin x + A\cos\omega s
]
其中，(c) 是阻尼系数，(A) 是驱动强度，(\omega) 是驱动频率。
操作步骤
-
将二阶方程转化为一阶常微分方程组：
[
\begin{cases}
x’ = y \
y’ = -cy - \sin x + A\cos\omega s
\end{cases}
] -
使用
ode45求解器求解微分方程。 - 绘制角速度随时间的变化曲线和相平面轨迹。
MATLAB代码
%drivenpendulum1.m
clc;clear;close;tic;
A=0.5;
y0=[0 0];
c=0.05;
w= 0.7;
ti=0; tf=300*pi;
M=round(2*tf);
tau = linspace(ti, tf, M);
drp= @(tau,y) [y(2);
-c*y(2) - sin(y(1)) + A*cos(w*tau)];
[t,Y]=ode45(drp,tau,y0);
x=Y(:,1);
yp=Y(:,2);
plot(t,yp,'linewidth',2);grid on;axis tight;
xlabel('t');ylabel('d\theta/dt');title(['A = ',num2str(A)])
xlim([.8*tf tf]);figure;
for i=floor(M/3):M
X(i)=x(i);
V(i)=yp(i);
end
plot(X,V,'linewidth',2);grid on;axis tight;%phase plane
title(['A = ', num2str(A)]);xlabel('\theta');ylabel('d\theta/dt');toc
4.3 双摆
双摆由两个摆锤连接而成，其运动具有很强的初始条件敏感性，难以得到解析解。双摆的运动可以用一组一阶常微分方程来描述。
操作步骤
- 定义双摆的运动方程。
-
使用
ode45求解器求解微分方程。 - 绘制摆角随时间的变化曲线。
MATLAB代码
%doublependulum1.m
clc;clear;close;
g = 9.81;
m = 1; L = 1; theta1 = 1*pi/4; theta2 = 2*pi/4;
t = linspace(0, 40, 800);
y0=[theta1, theta2, 0, 0];
dp = @(t,y) [6/(m*L^2) * (2*y(3)-3*cos(y(1)-y(2))*y(4)) / ...
(16-9*cos(y(1)-y(2))^2);
6/(m*L^2)*(8*y(4)-3*cos(y(1)-y(2))*y(3)) /(16-9*cos(y(1)-y(2))^2);
-1/2*m*L^2*((6/(m*L^2) * (2*y(3)-3*cos(y(1)-y(2))*y(4)) / ...
(16-9*cos(y(1)-y(2))^2))*( 6/(m*L^2)*(8*y(4)-3*cos(y(1)-y(2))*y(3)) /...
(16-9*cos(y(1)-y(2))^2))*sin(y(1)-y(2))+3*g/L*sin(y(1)));
-1/2*m*L^2*(-( 6/(m*L^2) * (2*y(3)-3*cos(y(1)-y(2))*y(4)) / ...
(16-9*cos(y(1)-y(2))^2))*(6/(m*L^2)*(8*y(4)-3*cos(y(1)-y(2))*y(3)) /...
(16-9*cos(y(1)-y(2))^2))*sin(y(1)-y(2))+g/L*sin(y(2)))];
[t,Y]=ode45(dp,t,y0);
z=[Y(:,1) Y(:,2)]*180/pi; % radians --> degrees
plot(t,z,'linewidth',2);grid on;axis tight;
legend('\theta1','\theta2');xlabel('t');ylabel('degrees');
总结
本文介绍了几种常见动力学系统的建模与分析方法，包括 FHN 神经元模型、混合槽与间歇反应器中的化学反应、四旋翼动力学以及摆问题。通过使用 MATLAB 的 ODE 求解器，我们可以对这些系统的动态行为进行数值模拟和分析。这些方法在工程、物理和生物学等领域都有广泛的应用。
表格总结
| 系统名称 | 主要方程 | 求解方法 | MATLAB 代码文件 |
|---|---|---|---|
| FHN 模型 | 一组常微分方程 |
ode45
|
FitzHughNagumo3.m
,
FitzHughNagumo4.m
|
| 混合槽 | (\frac{dM}{dt} = v_{in}C_{in} - v_{out}\frac{M}{V_0 + v_{in}t - v_{out}t}) |
ODE45
,
dsolve
|
odesolve_mix2.m
|
| 间歇反应器 | 一组反应速率方程 |
ode23
|
ode5.m
|
| 四旋翼 | 一组二阶常微分方程转化为一阶方程组 |
ode45
|
quadrotor.m
|
| 无阻尼单摆 | (\frac{d^2\theta}{dt^2} = -\omega^2\sin\theta) |
ode45
|
pendulum1.m
|
| 阻尼驱动摆 | (x’’ = -cx’ - \sin x + A\cos\omega s) |
ode45
|
drivenpendulum1.m
|
| 双摆 | 一组一阶常微分方程 |
ode45
|
doublependulum1.m
|
mermaid 流程图
graph TD
A[选择系统] --> B{FHN 模型}
A --> C{混合槽}
A --> D{间歇反应器}
A --> E{四旋翼}
A --> F{无阻尼单摆}
A --> G{阻尼驱动摆}
A --> H{双摆}
B --> B1[定义方程]
B1 --> B2[使用 ode45 求解]
B2 --> B3[绘制相平面和时域曲线]
C --> C1[定义方程]
C1 --> C2[使用 ODE45 或 dsolve 求解]
C2 --> C3[绘制糖质量曲线]
D --> D1[定义方程]
D1 --> D2[使用 ode23 求解]
D2 --> D3[优化反应速率常数]
E --> E1[转化为一阶方程组]
E1 --> E2[使用 ode45 求解]
E2 --> E3[绘制高度和角度曲线]
F --> F1[转化为一阶方程组]
F1 --> F2[使用 ode45 求解]
F2 --> F3[绘制摆角曲线]
G --> G1[转化为一阶方程组]
G1 --> G2[使用 ode45 求解]
G2 --> G3[绘制角速度和相平面曲线]
H --> H1[定义方程]
H1 --> H2[使用 ode45 求解]
H2 --> H3[绘制摆角曲线]
5. 各系统分析总结与对比
5.1 不同系统的求解方法对比
| 系统名称 | 方程类型 | 求解方法 | 特点 |
|---|---|---|---|
| FHN 模型 | 一阶常微分方程组 |
ode45
| 用于研究神经元放电行为的动态特性，通过改变激励电流观察极限环变化 |
| 混合槽 | 一阶常微分方程 |
ODE45
、
dsolve
| 可进行数值和解析求解，分析混合槽中糖质量随时间的变化 |
| 间歇反应器 | 一阶常微分方程组 |
ode23
| 用于优化化学反应的速率常数，以实现产品的最大收率 |
| 四旋翼 | 二阶常微分方程组转化为一阶 |
ode45
| 模拟四旋翼的运动，包括起飞频率的计算和运动轨迹的绘制 |
| 无阻尼单摆 | 二阶常微分方程转化为一阶 |
ode45
| 研究单摆的运动规律，不同摆角下周期的变化 |
| 阻尼驱动摆 | 二阶常微分方程转化为一阶 |
ode45
| 分析阻尼和驱动对摆运动的影响，观察混沌行为 |
| 双摆 | 一阶常微分方程组 |
ode45
| 模拟双摆的复杂运动，其运动对初始条件敏感 |
5.2 各系统的应用场景总结
- FHN 模型 ：在神经科学领域，用于模拟神经元的放电行为，帮助理解神经元的信号传递机制。
- 混合槽 ：在化工、食品等行业，用于模拟混合过程中物质的浓度变化，优化混合工艺。
- 间歇反应器 ：在化学工业中，用于优化化学反应的条件，提高产品的收率和质量。
- 四旋翼 ：在航空航天、物流等领域，用于模拟四旋翼飞行器的运动，进行飞行控制和路径规划。
- 无阻尼单摆 ：在物理学教学和研究中，用于演示简单的机械振动现象，验证物理理论。
- 阻尼驱动摆 ：在非线性动力学研究中，用于观察混沌现象，研究系统的稳定性和分岔特性。
- 双摆 ：在物理学和工程学中，用于研究复杂的机械运动，验证数值求解方法的有效性。
6. 实际应用中的注意事项
6.1 数值求解的精度问题
在使用
ode45
、
ode23
等求解器时，需要注意求解的精度。一般来说，求解器会自动调整步长以保证精度，但在某些情况下，可能需要手动调整步长或选择更合适的求解器。例如，对于一些刚性问题，
ode15s
可能比
ode45
更合适。
6.2 初始条件的选择
初始条件对系统的解有重要影响。在实际应用中，需要根据具体问题合理选择初始条件。例如，在四旋翼模型中，初始位置和速度的选择会影响飞行器的起飞和飞行轨迹。
6.3 模型的简化与实际情况的差异
在建立模型时，为了简化问题，往往会进行一些假设和简化。这些简化可能会导致模型与实际情况存在一定的差异。在应用模型结果时，需要考虑这些差异，并进行适当的修正。例如，在四旋翼模型中，忽略了空气阻力和其他外部干扰，实际飞行中这些因素可能会对飞行器的性能产生重要影响。
7. 拓展与展望
7.1 多系统耦合的研究
在实际应用中，很多问题涉及多个动力学系统的耦合。例如，在航空航天领域，飞行器的运动不仅受到自身动力学的影响，还受到大气环境、控制系统等多个因素的耦合作用。未来可以研究如何将不同的动力学系统进行耦合建模和求解，以更准确地模拟实际问题。
7.2 深度学习与动力学建模的结合
深度学习在数据处理和模式识别方面具有强大的能力。可以将深度学习方法与动力学建模相结合，利用深度学习算法从大量的实验数据中提取系统的动力学特征，提高模型的准确性和泛化能力。例如，在四旋翼飞行控制中，可以使用深度学习算法对飞行器的运动数据进行分析和预测，实现更智能的飞行控制。
7.3 实时仿真与控制
随着计算机技术的发展，实时仿真和控制变得越来越重要。在实际应用中，需要对动力学系统进行实时监测和控制。未来可以研究如何实现动力学系统的实时仿真和控制，提高系统的响应速度和稳定性。例如，在工业自动化生产中，实时仿真和控制可以帮助优化生产过程，提高生产效率和质量。
mermaid 流程图
graph LR
A[实际问题] --> B[建立动力学模型]
B --> C{选择求解方法}
C --> C1[数值求解]
C --> C2[解析求解]
C1 --> D[使用 ODE 求解器]
D --> E[调整参数和精度]
E --> F[分析结果]
F --> G{结果是否满意}
G -->|否| B[建立动力学模型]
G -->|是| H[应用结果]
C2 --> I[使用符号计算工具]
I --> J[化简和求解]
J --> F[分析结果]
总结
本文详细介绍了几种常见动力学系统的建模与分析方法，包括 FHN 神经元模型、混合槽与间歇反应器中的化学反应、四旋翼动力学以及摆问题。通过使用 MATLAB 的 ODE 求解器，我们可以对这些系统的动态行为进行数值模拟和分析。同时，我们还讨论了实际应用中的注意事项和未来的拓展方向。这些方法和技术在工程、物理和生物学等领域都有广泛的应用前景。在实际应用中，需要根据具体问题选择合适的模型和求解方法，并注意数值求解的精度、初始条件的选择以及模型与实际情况的差异。未来，多系统耦合的研究、深度学习与动力学建模的结合以及实时仿真与控制将是动力学建模领域的重要发展方向。
扫一扫
734
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
