Electronic – Can Op-Amp (LM358) works as an integrator without dual mode power supply

I'm not entirely sure how to answer the why question, but simulation reproduces the same result; first the big picture:

But if you zoom in there are tiny little overshoot peaks in there...

which are actually pretty close to those measured. TL081 model is from TI, which is actually a very basic one... yet it predicted the overshoot pretty well.

Alas while I'm not completely sure what actually causes them, thanks to discussion below with the OP, all fingers point to the slew rate of the TL081. If I lower the slopes of the square wave by 3 orders of magnitude (from 1ns to 1us), then the output spikes go away (in simulation), but I also get some DC bias on the output (once it stabilizes):

Yeah, the datasheet slew rate for TL081 is 13 V/us. So in 1ns (first simulation) it would only be able to swing its output by 13mV, not enough to react to 25x greater input swing. Basically the output gets briefly pulled (up or down) by the input via the feedback resistor before the opamp can react to correct that. Actually it almost works like a voltage doubler during those spikes. Another way to check this is to measure the time span of the spikes. In [the first] simulation those spikes last 240ns (=4.15MHz), which is comparable to the 3-4MHz bandwidth of the TL081 (I note that the old datasheet said 4MHz but the new one says 3MHz.)

The problem starts with the square wave.I'm not sure if it should go below 0 or if it's ok to let it vary between 0 and +5V for example.

Obviously, a DC offset in the input wave doesn't matter to the output – it won't pass through the transformer.

So you can use pretty much any periodic function. Since you said yours is a scientific question:

make yourself a plot of the spectrum of the periodic function you apply to the input (spectrum = Fourier transform!). You know that all periodic function have line spectra. The zero-mean square wave has lines in its spectrum spaced every 1/period – that's pretty handy, because that means that you can put your low pass somewhere betweeen 1/period and 2/period, which implies it can have pretty relaxed transition between pass- and stopband.

The triangle is just a square wave convoluted with itself. That's awesome – because convolution in time domain is multiplication in frequency domain, which means that the triangle wave's spectrum is just the squared spectrum of the square wave, so the same filter can be used. However, for practical reasons: It's a lot harder to produce a good triangle wave than a good square wave, and especially if you consider switching losses, fully switched semiconductors are usually much more efficient.

Also, from a practical point of view: your transformer already is an inductive device – and you might very well get away by adding a capacitor in series instead of having RC lowpasses on the input, which would give you a resonant circuit around the frequency you feed in, with (virtually) no power lost, solving the DC offset problems at the same time.