Duhamel积分和ode45

有老哥知道为什么用ode45和dumamel积分求解同一个问题结果却相差很远是什么问题呢？
%%ode
T0 = 1.0;
zeta = 0.05;
Tp = 0.5;
ag = 0.5 * 9.81;
tf = 5.0;
omega0 = 2*pi/T0;
m = 1;
k = m * omega0^2;
c = 2 * zeta * omega0 * m;
y0 = [0; 0];
tspan = [0 tf];
options = odeset('RelTol',1e-6,'AbsTol',1e-9); % 设置更严格的容差选项
[t, y] = ode45(@oscDynamics, tspan, y0, options);
figure;
plot(t, y(:,1), 'b-', 'LineWidth', 2);
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Response Using ode45 with Adjusted Tolerances');
grid on;
function dydt = oscDynamics(t, y)
zeta = 0.05;
T0 = 1.0;
omega0 =2*pi/T0;
m = 1;
k = m *omega0^2;
c = 2 * zeta *omega0 * m;
Tp= 0.5;
ag = 0.5 * 9.81;
dydt = zeros(2,1);
if t<= Tp
a = ag;
else
a = 0;
end
dydt(1) = y(2);
dydt(2) = (-c *y(2) - k * y(1) + m * a) / m;
end
%%
T0= 1.0;
zeta= 0.05;
Tp= 0.5;
ag= 0.5 * 9.81;
tf = 5.0;
omega0 = 2*pi/T0;
omega_d = omega0 * sqrt(1 - zeta^2);
m = 1;
k = m * omega0^2;
c = 2 * zeta * omega0 * m;
dt= 0.00001; %更小的时间步长以提高精度
t = 0:dt:tf;
ag_t = zeros(size(t));
ag_t(1:floor(Tp/dt)+ 1) = ag;
u = zeros(size(t));
exp_decay = exp(-zeta*omega0*dt);
sin_inc = sin(omega_d*dt);
cos_inc = cos(omega_d*dt);
s =0;
c =0;
for i = 2:length(t)
s_new = exp_decay * (cos_inc * s - sin_inc* c) + dt * ag_t(i);
c_new = exp_decay * (sin_inc * s + cos_inc* c);
s = s_new;
c = c_new;
u(i) = s / (omega_d * k);
end
figure;
plot(t,u, 'b-', 'LineWidth',2);
xlabel('Time (s)');
ylabel('Displacement (m)');
title('Numerical Response Using Enhanced Duhamel Integral');
gridon;