Solution of a Steady-State Heat Conduction Problem using the Galerkin Method  
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

Reference :
M. N. Özisik, Heat Conduction, Wiley, New York, 1980.

Approximate one-term trial Solution. We need to determine (Subscript[c, 1]).

T1[x_, y_] := (Subscript[c, 1]) (Subscript[ϕ, 1])[x, y]

(Subscript[ϕ, 1])[x_, y_] := (a^2 - x^2) (b^2 - y^2)

sol1 = 
      Solve[Integrate[
            Integrate[(D[T1[x, y], {x, 2}] + D[T1[x, y], {y, 2}] + 
                    1/k g) (Subscript[ϕ, 1])[x, y], {x, -a, a}], {y, -b, 
              b}] == 0, (Subscript[c, 1])]

plt1 = 
      Plot[T1[0, y] /. sol1 /. a -> 1 /. b -> 1 /. 
          g/k -> 1, {y, -1, 1}, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction ->  Identity]

plt2 = 
      Plot[T1[x, 0.5] /. sol1 /. a -> 1 /. b -> 1 /. 
          g/k -> 1, {x, -1, 1}, PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction ->  Identity]

Approximate two-term trial Solution. We need to determine (Subscript[c, 1])and (Subscript[c, 2]).

T2[x_, y_] := (Subscript[c, 1]) (Subscript[ϕ, 1])[x, y] + (Subscript[c, 2]) (Subscript[ϕ, 2])[x, y]

(Subscript[ϕ, 2])[x_, y_] := x^2 (a^2 - x^2) (b^2 - y^2)

sol2 = 
      Solve[{Integrate[
              Integrate[(D[T2[x, y], {x, 2}] + D[T2[x, y], {y, 2}] + 
                      1/k g) (Subscript[ϕ, 1])[x, y], {x, -a, 
                  a}], {y, -b, b}] == 0, 
          Integrate[
              Integrate[(D[T2[x, y], {x, 2}] + D[T2[x, y], {y, 2}] + 
                      1/k g) (Subscript[ϕ, 2])[x, y], {x, -a, 
                  a}], {y, -b, b}] == 0}, {(Subscript[c, 1]), (Subscript[c, 2])}]

plt3 = 
      Plot[T2[0, y] /. sol2 /. a -> 1 /. b -> 1 /. 
          g/k -> 1, {y, -1, 1}, PlotStyle -> RGBColor[0, 0, 1], 
        DisplayFunction ->  Identity]

plt4 = 
      Plot[T2[x, 0.5] /. sol2 /. a -> 1 /. b -> 1 /. 
          g/k -> 1, {x, -1, 1}, PlotStyle -> RGBColor[0, 0, 1], 
        DisplayFunction ->  Identity]

Exact Solution

β[n_] := (2 n + 1) Pi/2

Te[x_, y_] := 
      g/k (a^2 - x^2)/2  - 
        2 a^2 g/
            k (∑_(n = 0)^50)(-1)^
                  n/β[
                    n]^3 ((Cosh [β[n] (y/b)] Cos[β[n] (x/a)])/(Cosh [β[n] (b/a)] ))

plt5 = 
      Plot[Te[0, y] /. a -> 1 /. b -> 1 /. 
          g/k -> 1, {y, -1, 1}, PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction ->  Identity]

plt6 = 
      Plot[Te[x, 0.5] /. a -> 1 /. b -> 1 /. 
          g/k -> 1, {x, -1, 1}, PlotStyle -> RGBColor[0, 1, 0], 
        DisplayFunction ->  Identity]

The various solutions show good agreement.

T1[0, 0] /. sol1 /. a -> b // N

T2[0, 0] /. sol2 /. a -> b // N

Te[0, 0] /. a -> b // N // Simplify

Show[plt1, plt3, plt5, 
      DisplayFunction ->  $DisplayFunction]

Show[plt2, plt4, plt6, 
      DisplayFunction ->  $DisplayFunction]