3
$\begingroup$

I am attempting to find an efficient (less computationally demanding than solving with numeric integration) solution for the 1D heat equation:

$u_t - \alpha u_{xx} = 0$

IC: $ u(x, 0) = 0$

BC 1: $-u_x(0, t) = 1$

BC 2: $u(L, t) = 0$

I have attempted to solve this using separation of variables and have arrived at

$u(x, t) = L - x -\sum_{n=1}^{\infty}A_n e^{-\lambda_n^2\alpha t} \cos(\lambda_n x)$

with $ \lambda_n = \frac{(2n-1)\pi}{2L}, n = 1, 2, \dots $ and $ A_n = \frac{{(2n-1)(L-1)} \pi \left(-1\right)^{n} - 2}{L \lambda_{n}^{2}} $

Unfortunately, this series does not seem to converge. I can show my working to arrive at the above if it is required, but I feel certain that this problem has been solved before, so any pointers would be appreciated.

I include the solution steps I took below:

Let's say \begin{align} u(x, t) = w(x, t) + v(x) \end{align} where $w$ is the transient solution and $v$ is the steady state solution. Now, clearly $v(x) = c(L - x)$. This satisfies BC2 directly. BC1 leads to $-u_x(0, t) = -v_x = c = 1$, so \begin{align} v(x) = L-x \end{align}

For the transient solution, we go back to original problem, which transforms simply to $w_t = \alpha w_{xx}$.

IC: $u(x, 0) = 0 = w(x, 0) + v(x) = w(x, 0) + L - x$, so $w(x, 0) = x-L$

BC1: $-u_x(0, t) = -(w_x(0, t) + v_x(x)) = -w_x(0, t) + 1 = 1$, therefore $w_x(0, t)=0$.

BC2: $u(L, t) = w(L, t) + v(L) = w(L, t) = 0$.

We assume that $w(x, t) = X(x)T(t)$, so that $w_t=XT'$ and $w_{xx} = X''T$. The original problem therefore becomes $XT'=\alpha X''T$, which leads to \begin{align} \frac{X''}{X} = \frac{1}{\alpha}\frac{T'}{T} = -\lambda^2 \end{align} since these can only be equal if they are constant. Other signs or zero value for the constant leads only to trivial solutions.

We now have two separate ODEs: \begin{align} X'' + \lambda^2X = 0 \quad \text{and} \quad T + \lambda^2 \alpha T' = 0 \end{align} The general solutions to these are \begin{align} X = c_1 \cos (\lambda x) + c_2 \sin(\lambda x) \quad \text{and} \quad Y=c_3e^{-\lambda^2\alpha t} \end{align} Therefore \begin{align} w &= XT \\ &= c_3e^{-\lambda^2\alpha t}[ c_1 \cos(\lambda x) + c_2 \sin(\lambda x) ] \\ &= e^{-\lambda^2\alpha t}[A \cos(\lambda x) + B \sin(\lambda x)] \end{align} Applying the boundary conditions: \begin{align} BC1: w_x(0, t) = 0 \rightarrow -e^{-\lambda^2\alpha t}(A \lambda {\sin(\lambda 0)} + B \lambda \cos(\lambda 0)) = 0 \rightarrow B=0 \end{align} So we know that $w = Ae^{-\lambda^2\alpha t}\cos(\lambda x)$ \begin{align} BC2: w(L, t) = 0 \rightarrow w = Ae^{-\lambda^2\alpha t}\cos(\lambda L) = 0 \end{align} Obviously there are infinitely many $\lambda$s that satisfy this condition, which we indicate by writing \begin{align} \lambda_n = \frac{(2n-1)\pi}{2L}, n = 1, 2, \dots \end{align} Since all these values are solutions to the original problem, but do not in general satisfy the boundary conditions individually, we assume that some linear combination of the solutions is the actual solution, and write \begin{align} w = \sum_{n=1}^\infty A_ne^{-\lambda_n^2\alpha t}\cos(\lambda_nx) \end{align}

We solve for $A_n$ from the initial condition: \begin{align} w(x, 0) = x - L = \sum_{n=1}^\infty A_n\cos(\lambda_nx) \end{align}

If we multiply by $\cos(\lambda_m x)$ (note the different subscript) and integrate along the length of the rod (0 to $L$), we find \begin{align} \int_0^L \cos(\lambda_m x)x - L\cos(\lambda_mx) \mathrm{d}x = A_n \sum_{n=1}^\infty \int_0^L \cos(\lambda_nx) \cos(\lambda_m x) \mathrm{d}x \end{align} The integral on the right vanishes except when $n=m$, so we have \begin{align} \int_0^L \cos(\lambda_n x)x\mathrm{d}x - L\int_0^L\cos(\lambda_nx) \mathrm{d}x = A_n \sum_{n=1}^\infty \int_0^L \cos^2(\lambda_nx) \mathrm{d}x \end{align} Which evaluates to \begin{align} \frac{L \lambda_{n} \sin(L \lambda_{n}) + \cos(L \lambda_{n}) - 1}{\lambda_{n}^{2}}% -\frac{L^{2} \sin(L \lambda_{n})}{\lambda_{n}}% =% A_n\frac{2L\lambda_n + \sin(2L\lambda_n)}{4\lambda_n} \end{align} We can solve for $A_n$ as \begin{align} A_n = -\frac{4 \, {\left({\left(L^{2} \lambda_{n} - L \lambda_{n}\right)} \sin\left(L \lambda_{n}\right) - \cos\left(L \lambda_{n}\right) + 1\right)}}{2 \, L \lambda_{n}^{2} + \lambda_{n} \sin\left(2 \, L \lambda_{n}\right)} \end{align} Applying trigonometric identites yields \begin{align} A_n = -\frac{2 \, {\left({\left(L^{2} \lambda_{n} - L \lambda_{n}\right)} \sin\left(L \lambda_{n}\right) - \cos\left(L \lambda_{n}\right) + 1\right)}}{L \lambda_{n}^{2} + \lambda_{n} \sin\left(L \lambda_{n}\right) \cos\left(L \lambda_{n}\right)} \end{align} But we found $\lambda_n$ by requiring $\cos(L\lambda_n)=0$, so \begin{align} A_n = -\frac{2 \, {\left({\left(L^{2} \lambda_{n} - L \lambda_{n}\right)} \sin\left(L \lambda_{n}\right) + 1\right)}}{L \lambda_{n}^{2}} \end{align}

Now we substitute our solution for $\lambda_n$ and simplify to find \begin{align} A_n = \frac{{(2n-1)(L-1)} \pi \left(-1\right)^{n} - 2}{L \lambda_{n}^{2}} \end{align}

$\endgroup$
  • $\begingroup$ In the fifth equation from the bottom, the second term on the LHS is $$\frac{L \sin{\left (L \lambda_n\right)}}{\lambda_n}$$ $\endgroup$ – Ron Gordon May 25 '13 at 19:38
  • $\begingroup$ That's the integral, but remember that there is an L in front of it. $\endgroup$ – chthonicdaemon May 26 '13 at 7:15
  • $\begingroup$ $$\int_0^L \cos{\lambda_n x} \, dx = \frac{1}{\lambda_n} \sin{\lambda_n L}$$ $\endgroup$ – Ron Gordon May 26 '13 at 12:43
  • $\begingroup$ Indeed! When I fix that error and work through the rest of the problem, I get the same as your solution. Thank you very much for your help. $\endgroup$ – chthonicdaemon May 26 '13 at 14:13
2
$\begingroup$

I'm sorry, but I cannot tell where you went wrong; your answer makes no sense from a dimensional standpoint. As for the convergence of the series, I will answer below, but first let's actually solve the equation correctly.

EXACT SOLUTION OF THE EQUATION

In my humble opinion, I think these particular mixed BCs do not lend themselves easily to a simple solution by separation of variables (although the final solution is indeed of that form). I chose, rather, to solve this equation using a Laplace transform over $t$; that is

$$\hat{u}(x,s) = \int_0^{\infty} dt \, u(x,t) e^{-s t}$$

Then the PDE becomes on ODE on $x$:

$$\alpha \frac{d^2}{d x^2} \hat{u}(x,s) - s \hat{u}(x,s)=0$$

$$\left[\frac{d}{dx} \hat{u}(x,s)\right]_{x=0} = -\frac{1}{s}$$

$$\hat{u}(L,s) = 0$$

Note that the ODE already incorporates the initial condition.

The solution to this equation (the details of which solution I leave to the reader) is

$$\hat{u}(x,s) = \sqrt{\alpha} \frac{\sinh{\left[\frac{L-x}{\sqrt{\alpha}} \sqrt{s}\right]}}{s^{3/2} \cosh{\left[\frac{L}{\sqrt{\alpha}} \sqrt{s}\right]}}$$

Getting the solution to the original equation is then a matter of taking an inverse Laplace transform:

$$u(x,t) = \frac{1}{i 2 \pi} \int_{c-i \infty}^{c+i \infty} ds \, e^{s t} \hat{u}(x,s)$$

Although the integration contour appears to be complicated by the branch points in the LT expression, the ILT turns out to be simply equal to the sum of the residues of the poles of $\hat{u}(x,s)$. The reason this is so is because the functions of $\sqrt{s}$ are all even, so that there is no effect of crossing any branch cuts. The poles of $\hat{u}$ are at $s = -(n+1/2)^2 \pi^2 \alpha/L^2$ for integers $n$ and at $s=0$. Using the fact that, for a function of the form $f(s) = p(s)/q(s)$, the residue of a simple pole $s_0$ of $f$ is $p(s_0)/q'(s_0)$, I get for the ILT, after some algebra:

$$u(x,t) =L-x- 2 \frac{L}{\pi^2} \sum_{n=0}^{\infty} \frac{1}{(n+1/2)^2} \cos{\left[ \left ( n+\frac12 \right) \pi \frac{x}{L}\right]} e^{-\alpha (n+1/2)^2 \pi^2 t/L^2}$$

The series obviously converges, and you can perform some rudimentary checks to see that the initial and boundary conditions are satisfied, e.g., note that

$$\sum_{n=0}^{\infty} \frac{1}{(n+1/2)^2} = \frac{\pi^2}{2}$$

Also, you should check that the solution jibes with intuition at the steady state (i.e., as $t \to \infty$): here, $u(x,t) \sim L-x$, which makes physical sense (the temperature gradient at $x=0$ implies heat flow out of the left end of the bar).

RATE OF CONVERGENCE

The series above, while clearly convergent, does so painfully slow, especially for small $t$. Even for $n=50$, there is significant oscillation in the sum for $t=0$ which should ultimately converge to zero everywhere in $(0,L)$. Nevertheless, as $t$ increases, I do see increasingly fast convergence. The following is the behavior of the solution for $L=\alpha=1$, $n=10$ terms for times between $t=0.01$ and $t=5$:

enter image description here

$n=20$ is identical.

$\endgroup$
  • $\begingroup$ Thank you very much for this. Those temperature profiles look very much like what I would expect. I have expanded the original question to include my working, and I would appreciate if you had a look over it and point out where I went wrong. From what I can tell the problem came in at the calculation of $A_n$, as the rest of my solution seems to correspond with yours. $\endgroup$ – chthonicdaemon May 25 '13 at 19:29
  • $\begingroup$ @chthonicdaemon: my pleasure. I did point out a mistake just now in your derivation. The equation was dimensionally correct, so your answer needs to be. I should have thought to subtract out a transient term, this was a trick I learned in my Analysis class many moons ago. In a way, I am glad I didn't because I got to apply Laplace transforms which you will see is a very powerful and more general technique for solving PDE's over semi-infinite intervals. The transient term here comes from the pole at $s=0$. $\endgroup$ – Ron Gordon May 25 '13 at 19:46
  • $\begingroup$ One troubling element of this solution is that it shows energy going the wrong way. The physical situation is an imposed flux at x=0 and a fixed end point temperature, but the temperature goes below the initial value somehow for small $t$, even though there is no colder source available. I suppose this just shows that numerical solution is the right way to go here. $\endgroup$ – chthonicdaemon May 26 '13 at 12:03
  • $\begingroup$ @chthonicdaemon: forgive me, but I have no idea what you are saying. The condition imposed at $x=0$ is that the gradient is $-1$ for all $t>0$. As you can see from the above plot of the solution, that is certainly the case. (Note how the curves are parallel at $x=0$.) Any other physical features (e.g., "energy", "heat flux", etc.) may or may not obey your physical intuition, because they do not have to. What is a bit crazy about the analytical solution is how it behaves near $x=0$ for small $t$: the imposition of the BC there causes the series to oscillate badly and not converge well. $\endgroup$ – Ron Gordon May 26 '13 at 12:33
  • $\begingroup$ I'm sorry - I was evaluating the answer using code I had written for mine, and had not caught on that your n started at 0 while mine started at 1, so I wasn't adding in the first term, which made things look quite different. I have fixed it now and your answer is correct as far as I can see. $\endgroup$ – chthonicdaemon May 26 '13 at 13:47

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.