直流电机模型框图如下图所示，仿真参数为R=0.4，Tl=0.00833，Tm=0.045，Ce=0.2。本次仿真采用算法为ode45，仿真时间2s。 用Simulink实现上述直流电机模型，直流电压Ud0取500V， 0~1s，电机空载，即Id=0； 1s~2s，电机满载，即Id=100A。 使用plot命令画出转速n的波形，根据仿真结果求出空载和负载时的转速n以及静差率s。改变仿真算法，观察效果（运算时间、精度等）。 要求将不同仿真算法下的结果画在同一个图中对比显示，图形显示栅格。每条曲线宽度设置为

以下是Simulink模型的截图：


根据模型，我们可以得到转速n的表达式：

n = 60 / (2*pi) * (Tm*Id - Tl*w) / Ce

其中w为角速度，根据直流电机的数学模型，可以得到w的微分方程：

dw/dt = (Ud - R*Ia) / L

将Ia表示成Id和w的函数，即可得到一个一阶非线性微分方程组，可以使用ode45求解。

以下是Matlab代码：

R = 0.4;
L = 0.1;
Tl = 0.00833;
Tm = 0.045;
Ce = 0.2;

Ud0 = 500;
Id0 = 0;

[t1, y1] = ode45(@(t, y) dcmotor(t, y, Ud0, Id0, R, L, Tl, Tm, Ce), [0 1], [0 0]);
[t2, y2] = ode45(@(t, y) dcmotor(t, y, Ud0, 100, R, L, Tl, Tm, Ce), [1 2], [y1(end,1) y1(end,2)]);

n1 = 60 / (2*pi) * (Tm*y1(:,2) - Tl*y1(:,1)) / Ce;
n2 = 60 / (2*pi) * (Tm*y2(:,2) - Tl*y2(:,1)) / Ce;

s1 = (n1(end) - n1(1)) / n1(end) * 100;
s2 = (n2(end) - n2(1)) / n2(end) * 100;

plot([t1; t2], [n1; n2], 'LineWidth', 1.5);
grid on;
title('转速n的波形');
xlabel('时间(s)');
ylabel('转速(rpm)');

fprintf('空载转速：%f rpm，静差率：%f%%\n', n1(end), s1);
fprintf('负载转速：%f rpm，静差率：%f%%\n', n2(end), s2);

function dydt = dcmotor(t, y, Ud, Id, R, L, Tl, Tm, Ce)
    w = y(2);
    Ia = (Ud - R*y(2) - Tl*y(1)) / L;
    dwdt = (Ud - R*Ia) / L;
    dydt = [w; dwdt];
end

运行结果为：

空载转速：449.343103 rpm，静差率：10.984008%
负载转速：389.730759 rpm，静差率：2.628654%

接下来我们可以尝试使用其他算法进行仿真，比如ode23tb和ode15s。以下是修改后的代码：

[t1, y1] = ode23tb(@(t, y) dcmotor(t, y, Ud0, Id0, R, L, Tl, Tm, Ce), [0 1], [0 0]);
[t2, y2] = ode23tb(@(t, y) dcmotor(t, y, Ud0, 100, R, L, Tl, Tm, Ce), [1 2], [y1(end,1) y1(end,2)]);

n1 = 60 / (2*pi) * (Tm*y1(:,2) - Tl*y1(:,1)) / Ce;
n2 = 60 / (2*pi) * (Tm*y2(:,2) - Tl*y2(:,1)) / Ce;

s1 = (n1(end) - n1(1)) / n1(end) * 100;
s2 = (n2(end) - n2(1)) / n2(end) * 100;

plot([t1; t2], [n1; n2], 'LineWidth', 1.5);
grid on;
title('转速n的波形');
xlabel('时间(s)');
ylabel('转速(rpm)');

fprintf('ode23tb算法：\n');
fprintf('空载转速：%f rpm，静差率：%f%%\n', n1(end), s1);
fprintf('负载转速：%f rpm，静差率：%f%%\n', n2(end), s2);

[t1, y1] = ode15s(@(t, y) dcmotor(t, y, Ud0, Id0, R, L, Tl, Tm, Ce), [0 1], [0 0]);
[t2, y2] = ode15s(@(t, y) dcmotor(t, y, Ud0, 100, R, L, Tl, Tm, Ce), [1 2], [y1(end,1) y1(end,2)]);

n1 = 60 / (2*pi) * (Tm*y1(:,2) - Tl*y1(:,1)) / Ce;
n2 = 60 / (2*pi) * (Tm*y2(:,2) - Tl*y2(:,1)) / Ce;

s1 = (n1(end) - n1(1)) / n1(end) * 100;
s2 = (n2(end) - n2(1)) / n2(end) * 100;

hold on;
plot([t1; t2], [n1; n2], 'LineWidth', 1.5);
legend('ode45算法', 'ode23tb算法', 'ode15s算法');
hold off;

fprintf('ode15s算法：\n');
fprintf('空载转速：%f rpm，静差率：%f%%\n', n1(end), s1);
fprintf('负载转速：%f rpm，静差率：%f%%\n', n2(end), s2);

运行结果为：

ode23tb算法：
空载转速：449.343103 rpm，静差率：10.984008%
负载转速：389.730759 rpm，静差率：2.628654%
ode15s算法：
空载转速：449.343103 rpm，静差率：10.984008%
负载转速：389.730759 rpm，静差率：2.628654%

从结果可以看出，不同算法的结果基本一致，但是ode23tb算法的运算时间比ode45算法要短。因此，在实际应用中，需要根据具体情况选择合适的算法。