I'll try to answer to your questions in their order of appearance below.
- How is the variation $\delta y$ different from the differential $dy$? one of the books i tried to use to understand this says the former is "virtual" while the later is "real". Now I am an engineer so this might sound like a dumb thing to say among mathematicians, but they both sound like they mean the same thing. Maybe there is a formality here that I am not understanding. Are differentials a special case of variations ?
Technically speaking, differentials are variations of function(al)s defined of domains on $\Bbb R$ or $\Bbb C$ (or more generally finite dimensional vector spaces), thus they are special cases of the variation concept. In this respect, perhaps it would also be useful to have a look at this Q&A on functional derivatives: the variation $\delta T$ of the argument of your functional is not a simple variable but a function. For example, in our case,
$$
\delta T=\varepsilon \theta \quad\varepsilon\in]0,1]
$$
where $\theta$ is a relatively small (respect to the sought for stationary $T$ value) but admissible temperature field increment, meaning with this that $I(T+\theta)$ is perfectly defined. "Virtuality" is an attribute given to variations that arises from the application of the principle of virtual work, itself an application of the variational principle of least action: however, they are simply variations.
- Why wasn't the Euler-Lagrange used here ? Is it because we've now discretized the integral and it is now a summation ? Or is there another reason ?
The Euler Lagrange-Equation is not used here since it is useless in this context: the basic aim of the method of finite elements (and likewise of the finite difference method, despite their large differences) is to approximate a vector field with a piecewise constant vector function (a simple function using the terminology of integration theory) obtained by solving a linear algebraic system, hopefully sufficiently small to be tractable. In the case under analysis, you proceed by dividing your domain in a (finite) number $n$ of regions (elements) and evaluate the given functional $I(T)$ in each of them by doing the following steps in this exact order
- Approximate the temperature field $T$ inside each region by a piecewise constant function, and then
- Minimize the local functional, now transformed in a simpler linear algebraic relation by calculating its functional derivative and impose its vanishing.
If you apply the functional derivative before choosing an approximation for the temperature field, you get exactly the Euler-Lagrange equation of the functional $I(T)$ but, as you can check, this is simply the local formulation of the former problem you are trying to solve, i.e. the stationary equation for the heat conduction
$$
\frac{\partial}{\partial x}\left(k_x\frac{\partial T}{\partial x}\right)
+ \frac{\partial}{\partial y}\left(k_y\frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial z}\left(k_z\frac{\partial T}{\partial z}\right) = G\quad \quad (x,y,z)\in \Omega_e,\; e=1, \dots, M,
$$
therefore you haven't done any step forward towards the resolution of your problem (but possibly a backward step).
