In even space-time dimensions the multi-loop Feynman integrals are integrals
of rational function in projective space. By using an algorithm that extends
the Griffiths–Dwork reduction for the case of projective hypersurfaces with
singularities, we derive Fuchsian linear differential equations, the
Picard–Fuchs equations, with respect to kinematic parameters for a large class
of massive multi-loop Feynman integrals. With this approach we obtain the
differential operator for Feynman integrals to high multiplicities and high
loop orders. Using recent factorisation algorithms we give the minimal order
differential operator in most of the cases studied in this paper. Amongst our
results are that the order of Picard–Fuchs operator for the generic massive
two-point $n-1$-loop sunset integral in two-dimensions is
$2^{n}-\binom{n+1}{\left\lfloor \frac{n+1}{2}\right\rfloor }$ supporting the
conjecture that the sunset Feynman integrals are relative periods of
Calabi–Yau of dimensions $n-2$. We have checked this explicitly till six
loops. As well, we obtain a particular Picard–Fuchs operator of order 11 for
the massive five-point tardigrade non-planar two-loop integral in four
dimensions for generic mass and kinematic configurations, suggesting that it
arises from $K3$ surface with Picard number 11. We determine as well
Picard–Fuchs operators of two-loop graphs with various multiplicities in four
dimensions, finding Fuchsian differential operators with either Liouvillian or
elliptic solutions.