2D impulse physics engine broken polygon restitution

Been trying to make a simple 2D physics engine in C#. Followed a tutorial by "Two-Bit Coding" on youtube but noticed that squares, with two contact points would not bounce at all. Tracked that problem down to the impulses simply being too small when resolving the collisions.

Got some bounciness by multiplying the impulse by a value. This however is not a solution, since if the collision occurs very close to the body-center (on the middle an edge of a square) the impulse would be massive, causing the bounce to break the laws of energy conservation.

Tried another formula to resolve the collision, this time from Wikipedia.The formula is for 3D but I assumed it would work for 2D. Formula from wikipedia

This is the C# implementation:

        public void Resolve_collision_wiki(in Physics_manifold contact)
        {
            Physics_body body_a = contact.body_a;
            Physics_body body_b = contact.body_b;
            Physics_vector normal = contact.normal;

            int contact_count = contact.contact_count;

            float e = 1;

            contact_list[0] = contact.contact_point_1;
            contact_list[1] = contact.contact_point_2;

            for (int i = 0; i < contact_count; i++)
            {
                Physics_vector ra = contact_list[i] - body_a.Position;
                Physics_vector rb = contact_list[i] - body_b.Position;

                center_a_contact_list[i] = ra;
                center_b_contact_list[i] = rb;


                Physics_vector relative_velocity = body_b.Linear_velocity - body_a.Linear_velocity + Physics_math.Cross_product(body_b.Angular_velocity, rb) - Physics_math.Cross_product(body_a.Angular_velocity, ra);

                //() first, then the dot product?
                float numerator = Physics_math.Dot_product(-(1 + e) * relative_velocity, normal);

                if (numerator < 0)
                    return;

                //https://www.cs.cmu.edu/~baraff/sigcourse/notesd2.pdf#page=16

                numerator /= body_a.inverse_mass + body_b.inverse_mass + Physics_math.Dot_product(Physics_math.Cross_product(body_a.inverse_inertia * Physics_math.Cross_product(ra, normal), ra) + Physics_math.Cross_product(body_b.inverse_inertia * Physics_math.Cross_product(rb, normal), rb), normal);
                numerator = MathF.Max(numerator, 0);

                j_list[i] = numerator;
            }


            //Applying in the first loop would make the second impulse be a bit off
            for (int i = 0; i < contact_count; i++)
            {
                body_a.Linear_velocity -= j_list[i] * body_a.inverse_mass * normal;
                body_a.Angular_velocity -=  Physics_math.Cross_product(center_a_contact_list[i], normal * j_list[i]) * body_a.inverse_inertia;

                body_b.Linear_velocity += j_list[i] * body_b.inverse_mass * normal;
                body_b.Angular_velocity +=  Physics_math.Cross_product(center_b_contact_list[i], normal * j_list[i]) * body_b.inverse_inertia;
            }
        }

The dot product and 2D cross product implementation:

        public static float Dot_product(Physics_vector a, Physics_vector b)
        {
            return a.x * b.x + a.y * b.y;
        }

        public static float Cross_product(Physics_vector a, Physics_vector b)
        {
            return a.x * b.y - a.y * b.x;
        }

        public static Physics_vector Cross_product(float a, Physics_vector v)
        {
            return new Physics_vector(-a * v.y, a * v.x);
        }

But the bounce height still changes dramatically when I move to contact points closer to each other, and gets lower when the contact points are farther apart. Furthermore, an e value of 1 won't give a perfectly elastic collision for squares like it does for circles.

Each contact is a separate velocity constraint that needs solving. The solution of these contact constraints must be simultaneous - the effect of the solution of one constraint (for contact constraints, the impulse applied at the contact point) on the other constraints must be accounted for.

To accurately solve each constraint we must thus solve a system of equations for all the active constraints at a given point in time. However, in practice solving this system of equations (which could contain thousands of equations) is very computationally expensive so realtime physics engines opt for an iterative approximation.

This iterative constraint solver, which you are trying to implement, solves each constraint individually over a number of iterations. To correctly approximate the simultaneous solution of all the constraints, the iterative solver must take into account the partial solution of the constraints it has already solved. Otherwise, the effect of each constraint on the others is overlooked.

Your implementation solves each constraint independently, and disregards the effect of the second constraint on the first (and vice versa). This is fine for circles since typical circle collision routines (presumably yours included) only generate a single contact to solve. However, this breaks down and gives an inaccurate solution when a collision generates multiple contact points. This is precisely why you have to apply the impulse in the first loop itself, contrary to your comment (the second impulse being "off" is it being affected by the result of the first, which is what we are looking for).

For even more stability, you could choose to iteratively solve all the contact constraints multiple times (taking into account the previous solutions) - the added iterations will help contact constraints across different pairs of bodies factor in their effects on each other. For instance, this will help noticeably reduce the jitter that can be seen when one body is sandwiched between two other bodies under the influence of gravity.