Peak gain on Bode plot not equals h infinity norm?

Hello I am using norm(H,inf) to calculate the h-infinity norm for a discrete-time system. But when I have a code like this:

load num, den;

H = tf(num,den,60);

bode(H)

[ninf,fpeak] = norm(H,inf);

x = 20*log10(ninf);

I found that x=38.2 while the peak gain on bode plot is 6.8. Why are they different? Is it because norm(X,inf) is not accurate? Btw, I cannot use hinfnorm in my case because it is a discrete-time system.

I have attached the code and .mat files.

Thanks in advance.

[Pascal Gahinet]

This is a spectacular illustration of the numerical pitfalls of high-order transfer functions discussed here:

https://www.mathworks.com/help/control/examples/using-the-right-model-representation.html

You can see there is a problem by plotting the singular values of the TF, ZPK, and SS representations of this model:

sigma(H,zpk(H),ss(H))

The response is all over the place for w<0.2 rad/s. PRESCALE confirms that the state-space representation has very high sensitivity below this frequency:

prescale(ss(H))

The system H is a transfer function of order 47 with an interesting pole/zero distribution:

pzmap(H)

Numerically this pattern is consistent with a pole/zero of high multiplicity, see also

https://www.mathworks.com/help/control/examples/sensitivity-of-multiple-roots.html

Such configuration is extremely sensitive to numerical errors, which essentially means that double precision is not enough to capture the system behavior at low frequency, and we really don't know what its frequency response might be, let alone its peak gain.

If possible, use the ZPK or SS representations to incrementally build H and accurately keep track of its dynamics, and stay away from the TF representation.