The governing differential equation is

$$\nabla^2 T=0 \tag A$$

The boundary conditions for this problem are as foll0ws:

$$T(0,y,z)=T_{hi} \tag {1A}$$

$$T(L,y,z) = T(x,0,z) = T(x,l,z) = T(x,y,0)= T(x,y,\mu) = 0 \tag {1B}$$

I need to solve for the distribution $T(x,y,z)$. I include my attempt here.

Solution attempt

Homogeneous Dirichlet type B.C. on $y$ and $z$ faces allow us to write the following form of preliminary temperature distribution:

$$T(x,y,z) = \sum_{n,m=1}^{\infty}\sin\bigg(\frac{n\pi y}{l}\bigg)\sin\bigg(\frac{m\pi z}{\mu}\bigg)T_{nm}(x) \tag {1C}$$

We substitute $(1\mathrm{C})$ in $\mathrm{A}$ and apply the orthogonality properties of $\sin\bigg(\frac{n\pi y}{l}\bigg)$ and $\sin\bigg(\frac{m\pi z}{\mu}\bigg)$ to obtain the following:

$$\frac{\mathrm{d^2}T_{nm}(x)}{\mathrm{d}x^2} - \gamma_{nm}^2 T_{nm}(x) = 0 \tag {1D}$$ where $\gamma_{nm} = \sqrt{(\frac{n\pi}{l})^2+(\frac{m\pi}{\mu})^2}$

The general solution of $(1\mathrm{D})$ is of the form:

$$T_{nm}(x) = A_{nm}e^{\gamma x} + B_{nm}e^{-\gamma x} \tag {1E}$$

$A_{nm},B_{nm}$ are the unknown Fourier coefficients that need to be determined. Applying $T(L,y,z) = 0$ on $(1\mathrm{D})$ and using $(1\mathrm{E})$, we arrive at:

$$A_{nm}e^{\gamma L} + B_{nm}e^{-\gamma L} = 0$$

$$\bf{B_{nm} = - A_{nm}e^{2\gamma L} \tag {1F}}$$

$$T_{nm}(x) = A_{nm}(e^{\gamma x} - e^{2\gamma L - \gamma x}) \tag {1G}$$

Using the B.C. $(1\mathrm{A})$:

$$T_{hi} = \sum_{n,m=1}^{\infty}\sin\bigg(\frac{n\pi y}{l}\bigg)\sin\bigg(\frac{m\pi z}{\mu}\bigg)A_{nm}(1 - e^{2\gamma L}) \tag {1H}$$

Multiplying both sides of $(1\mathrm{H})$ by $\int_0^l\sin\bigg(\frac{k\pi y}{l}\bigg)\mathrm{d}y$ and $\int_0^{\mu}\sin\bigg(\frac{j\pi z}{\mu}\bigg)\mathrm{d}z$ and using the principle of orthogonality, we arrive at:

$$T_{hi}\int_0^l \int_0^\mu \sin\bigg(\frac{k\pi y}{l}\bigg)\sin\bigg(\frac{j\pi z}{\mu}\bigg) \mathrm{d}y\mathrm{d}z = \frac{l\mu}{4} A_{kj} (1 - e^{2\gamma L}) \tag {1I}$$

Solving the definite integrals involved , we arrive at (for any arbitrary integer $n$ and $m$):

$$\bf{A_{nm} = \frac{4 T_{hi}}{nm\pi^2(1-e^{2\gamma L}) }(1-\cos(n\pi))(1-\cos(m\pi)) \tag {1J}}$$

$$T(x,y,z) = \sum_{n,m=1}^{\infty}(A_{nm}e^{\gamma x} + B_{nm}e^{-\gamma x})\sin\bigg(\frac{n\pi y}{l}\bigg)\sin\bigg(\frac{m\pi z}{\mu}\bigg) \tag {1K}$$

On coding this solution in MATLAB. and substituting $x=0$, I find that the answer is not equal to $T_{hi}$. For parameter values $L=0.5,l=0.5,\mu=0.05,T_{hi}=50$. When i evaluate for $x=0,y=l/2,z=\mu /2$, the answer should be $T=50$, but it evaluates to $96$ using the first four terms in the series. The series surely converges.

Is there something wrong, in the way, I am doing the problem ?Any help is greatly appreciated.

An observation

When I take $z=\mu /2, y=l/2$ and $x=0$, and consider only odd values for $n$ and $m$, the solution can be written as:

$$T\bigg(0,\frac{l}{2},\frac{\mu}{2}\bigg) = \frac{16 T_{hi}}{\pi^2}\underbrace{\bigg[1 + \frac{1}{9} + \frac{1}{25} + \frac{1}{49} + ........\bigg]}_{\pi^2 /8} = 2T_{hi} \tag {1L}$$

So, is there something wrong with my analytical solution ?


Are you sure you typed it in correctly? I checked your solution against mine and they're the same. I'll leave my answer below in case it helps anyways.

Often times I find it better to start from the beginning when I'm stuck to remove doubt of errors made in assumptions. \begin{cases} \begin{align} \Delta u(x,y,z) &= 0 \\ u(0,y,z) &= u_0 \\ u\vert_{\partial\Omega\backslash\{x=0\}} &= 0. \end{align} \end{cases} You know about separation of variables already, so we may take the eigenfunctions in the $y$ and $z$ directions as $\sin \frac{n\pi y}{\ell}$ and $\sin \frac{m\pi z}{\mu}$. The equation in $x$ is known as $a_{mn}e^{\gamma_{mn}x} + b_{mn}e^{-\gamma_{mn}x}.$ Thus $u$ is the resulting infinite linear combination $$ u(x,y,z) = \sum_{m,n \geq 1}\left(a_{mn}e^{\gamma_{mn}x} + b_{mn}e^{-\gamma_{mn}x}\right)\sin\frac{n\pi y}{\ell}\sin\frac{m\pi z}{\mu}. $$ The boundary conditions on $y$ and $z$ are already satisfied so we look towards those for $x$. \begin{cases} \begin{align} u_0 &= \sum_{m,n \geq 1}(a_{mn} + b_{mn})\sin\frac{n\pi y}{\ell}\sin\frac{m\pi z}{\mu} \\ 0 &= \sum_{m,n \geq 1} \left(a_{mn}e^{\gamma_{mn}L} + b_{mn}e^{-\gamma_{mn}L}\right)\sin\frac{n\pi y}{\ell}\sin\frac{m\pi z}{\mu} \end{align} \end{cases} Note that using orthogonality is not multiplying each side by the respective integrals; it is multiplying each side by the eigenfunctions with dummy indices and then integrating. Continuing with this in mind, \begin{cases} \begin{align} \int_0^\ell \int_0^\mu u_0 \sin\frac{n\pi y}{\ell} \sin\frac{m\pi z}{\mu} \, \mathrm{d}z \,\mathrm{d}y &= \frac{\mu\ell}{4}(a_{mn} + b_{mn}) \\ 0 &= a_{mn}e^{\gamma_{mn}L} + b_{mn}e^{-\gamma_{mn}L} \end{align} \end{cases} The first integral is easily found and we have the system of equations $$ \begin{pmatrix} 1 & 1 \\ e^{\gamma_{mn}L} & e^{-\gamma_{mn}L} \end{pmatrix} \begin{pmatrix} a_{mn} \\ b_{mn} \end{pmatrix} = \begin{pmatrix} \frac{4\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{\pi^2mn}u_0 \\ 0 \end{pmatrix}. $$ Inverting to solve for $a_{mn}$ and $b_{mn}$ then, \begin{align} \begin{pmatrix} a_{mn} \\ b_{mn} \end{pmatrix} &= -\frac{\operatorname{csch}(\gamma_{mn}L)}{2} \begin{pmatrix} e^{-\gamma_{mn}L} & -1 \\ -e^{\gamma_{mn}L} & 1 \end{pmatrix} \begin{pmatrix} \frac{4\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{\pi^2mn} \\ 0 \end{pmatrix} \\ &= -\frac{2u_0}{\pi^2}\frac{\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{\sinh(\gamma_{mn}L) mn} \begin{pmatrix} e^{-\gamma_{mn}L} \\ -e^{\gamma_{mn}L} \end{pmatrix}. \end{align} Thus $u$ is written \begin{equation} u(x,y,z) = -\frac{4u_0}{\pi^2} \sum_{m,n \geq 1} \Gamma_{mn} \frac{\sinh\big(\gamma_{mn}(x-L)\big)}{\sinh(\gamma_{mn}L)} \sin\frac{n \pi y}{\ell} \sin \frac{m\pi z}{\mu} \end{equation} with $$ \Gamma_{mn} = \frac{\big((-1)^m - 1\big) \big((-1)^n - 1\big)}{mn}. $$

$\hskip 1 in$ Laplace Eq sol

Mathematica code - it should be copy-pastable - enjoy.

\[Gamma]mn = Sqrt[((n \[Pi])/\[ScriptL])^2 + ((m \[Pi])/\[Mu])^2];
(* u0 not specified yet *)
L = 1;
\[ScriptL] = 1;
\[Mu] = 1;
(* Keep q low (~10-15 or so) as the number of terms grows as q^2 *)
(* i.e. higher q gives more accuracy but takes longer to run in effect *)
q = 5;

(*The coefficients amn and bmn {amn,bmn}={a,b}/.Solve[{a+b\[Equal]4/\
\[Pi]^2(((-1)^m-1)((-1)^n-1))/(m n)u0,a Exp[\[Gamma]mn L]+b Exp[-\
\[Gamma]mn L]\[Equal]0},{a,b}][[1]];*)

u[x_, y_, 
   z_] = (-4 u0)/\[Pi]^2 Sum[
    Sum[(((-1)^m - 1) ((-1)^n - 1))/(m n) Sinh[\[Gamma]mn (x - L)]/
       Sinh[\[Gamma]mn L] Sin[n \[Pi] y/\[ScriptL]] Sin[
       m \[Pi] z/\[Mu]], {n, 1, q}], {m, 1, q}];

Plot3D[u[0, y, z]/u0, {y, 0, \[ScriptL]}, {z, 0, \[Mu]}, 
 PlotLabel -> 
  "Laplace's Equation Solution\nat x=0 with\n\[ScriptL]=1, \[Mu]=1, \
25 terms", 
 AxesLabel -> {"y", "z", "\!\(\*FractionBox[\(u\), \(u0\)]\)"}, 
 Boxed -> False]

(* Something extra - it's computationally expensive (3d), at least for my \
computer,  so let it run for a bit *)
 ContourPlot[u[x, y, z]/u0, {x, 0, L}, {y, 0, \[ScriptL]}, 
  ColorFunction -> "DarkRainbow"], {z, 0, \[Mu]}]

(* alt+.to quit if it takes too long & reduce q *)
 u[x, y, z]/u0, {x, 0, L}, {y, 0, \[ScriptL]}, {z, 0, \[Mu]}, 
 Mesh -> False, ColorFunction -> "Rainbow", 
 PlotLabel -> "Contours of Temperature", AxesLabel -> {"x", "y", "z"}]
  • $\begingroup$ Thanks for such a detailed and well-explained answer. It seems the answer i arrived at is same with what you get . Maybe I did some mistake while writing the MATLAB script. Saying that would it be possible for you to post your Mathematica script here ? Also, your note about multiplying with respective eigen functions was helpful. $\endgroup$ – Indrasis Mitra Mar 24 at 5:46
  • $\begingroup$ Sure, glad I could help. I already deleted the code but I could post just the solution and the plot if you'd like. $\endgroup$ – AEngineer Mar 24 at 5:58
  • $\begingroup$ That would be really helpful. Appreciate it. $\endgroup$ – Indrasis Mitra Mar 24 at 6:01
  • $\begingroup$ Thanks for the Mathematica script. Some really cool graphics included. I was hoping you could take a look at this problem which I have been trying recently, especially its SP3 part. Any guidance will be really helpful. $\endgroup$ – Indrasis Mitra Mar 24 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.