Resonance frequency from bode and damp do not agree. Matlab bug or user fault?

I am working on system identification and I have spotted a consistent divergence between the resonance frequency of the system obtained using bode(sys) and peak picking and the one obtained using the natural frequency and damping ratio from damp(sys). Below follows an example for a simulated mass, spring, damper system.

 clear all,close all,clc

 Mo = 146; % (kg) 1/4 mass of the chassis
 m = 94; % (kg)
 kMo = 5000; % (N/m)
 kmo = 20000; % (N/m)
 cM = 2500; % (N*s)/m

 ex = 1;
rng(ex)  % change seed number
M = randi([146 781]); % change mass randomly in the interval of the expected values form the manufaturer
ran = randi([-10 10]); % random number
kM = kMo + ran; % Randomly change the stifness of the suspension
km = kmo + ran*10; % Randomly change the stifness of the tire
% Transfer functions Laplace
%----------------------------------------
% Characteristic polynomial (denominator)
den=    [m*M/(m*M),    (cM*(m+M))/(m*M),   (kM*m+M*(km+kM))/(m*M),   (cM*km)/(m*M),    (km*kM)/(m*M)];
% Numerator of transfer function X(s)/R(s)
num(1,:) = [0, (km*cM)/(m*M), (km*kM)/(m*M), 0, 0];
% Numerator of transfer function X(s)/F(s)
num(2,:) = [m/(m*M), 0, km/(m*M), 0,  0];
% Define the transfer functions
G12=tf(num(1,:),den);
G22=tf(num(2,:),den);
bd = figure;
[mag,phase,w] = bode(G12);
nominal = plot(w/(2*pi),20*log10(squeeze(mag)),'b');
title('System bode')
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid on 
hold all
[Wn,z]=damp(G12);
wn=Wn(2)./(2*pi);
Wr = Wn(2).*sqrt(1-2*(z(2).^2)); %Resonance frequency (rad/s)
wr = Wr/(2*pi); % Resonance frequency (Hz)

If now you pick the peak at the magnitude of the bode plot you may see that the respective frequency is 1.464 Hz (resonance frequency), while the one obtained via damp.m is 0.5186 Hz (wr).

Is there something that I am doing wrong? Is this to be expected for some reason? Thank you in advance for your help.

Answer 1

[Arkadiy Turevskiy]

The formula for computing resonant frequency as a function of a natural frequency that you used is only valid for 2nd order system with no zeros. You have 4th order transfer function with 3 zeros. The formula is not valid anymore.

[Kenblack]

Thank you for your answer. Could you please provide a reference on this?

[Arkadiy Turevskiy]

All the derivations of the formulas you listed assume 2nd order system with no zeros. For example, take a look at Ogata, Modern Control Engineering. I have 3rd edition, there it is on pp. 482-483.