I have Scala code that recreates the Crank-Nicolson solutions for the diffusion equations, and matches 'Excel for Scientists and Engineers' (Joe Billo, Wiley).

However, I would like to be able to replicate the results from H. Asan's 'Numerical computation of time lags and decrement factors' (including the decrement factor with very thin materials) and also the multi-layer positioning results from the BS2013 paper 'Optimizing insulation-thermal mass layer distribution from maximum time lag and minimum decrement factor point of view'.

Both papers suggest that Crank-Nicolson is at the heart of the approach, however there is a complexity with the treatment of energy flow at the outer and inner surface, and with different material layers, that are not covered in the simple diffusion treatments that I found.

How do I set up the matrix - or how do I connect a series of such matrices together?

I had thought that perhaps a simple treatment would give a given wall thickness a behavior like a delay factor and an amplitude scaling - but it does not seem to be possible to do this and sum the phase delays and multiply the decrements, because the results above suggest that the ordering is not transitive.

How can I proceed - given that I'm a lot better at coding than solving PDEs.


1 Answer 1


The heat equation for non-constant coefficients, which is more or less what you're doing here, takes the form, $$ \frac{\partial\psi}{\partial t}=\nabla\cdot\left(\kappa\nabla\psi\right) $$ We can simplify by using a single dimension: $$ \frac{\partial\psi}{\partial t}=\frac{\partial}{\partial x}\left(\kappa\frac{\partial\psi}{\partial x}\right) $$ Rather than applying the derivative operator to the term in the parenthesis (to get a diffusive term & an advective term), we discretize here (though I heartily encourage you to work through the math to get this next equation) to obtain \begin{align} \frac{\psi^{n+1}-\psi^n}{dt}&=\frac1{2dx^2}\left(\kappa_{i+1/2}\left[\psi_{i+1}^{n+1}-\psi_i^{n+1}\right]-\kappa_{i-1/2}\left[\psi^{n+1}_i-\psi^{n+1}_{i-1}\right]\right)\\ &+\frac1{2dx^2}\left(\kappa_{i+1/2}\left[\psi_{i+1}^{n}-\psi_i^{n}\right]-\kappa_{i-1/2}\left[\psi^{n}_i-\psi^{n}_{i-1}\right]\right) \end{align} where $$ \kappa_{i+1/2}=\frac{1}{2}\left(\kappa_{i+1}+\kappa_i\right) $$ You can solve for the $\psi^{n+1}$ terms, $$ -\delta\kappa_{i+1/2}\psi^{n+1}_{i+1}+ \left(1+\delta\left(\kappa_{i+1/2}+\kappa_{i-1/2}\right)\right)\psi^{n+1}_i-\delta\kappa_{i-1/2}\psi^{n+1}_{i-1}=RHS $$ where $\delta=dt/2dx^2$ and $RHS$ contains the $\psi^n$ terms; this is solved as the normal constant-coefficient case. However, you would need to generate the $\kappa(x)\to\kappa_i$ mapping, just as you do for $\psi(x)\to\psi_i$.

You can test this by letting \begin{align} \kappa(x)&=1+\left(\frac x5\right)^2 \\ \psi(x,0)&=\begin{cases}15 & x<0 \\ 20 & x=0 \\ 25 & x>0\end{cases} \end{align} As you let $\psi(x,t)$ evolve, you'll find that the solution tends to the stationary form $$ \psi(x,t)=20\left(1+\frac1\pi\tan^{-1}\left(\frac x5\right)\right) $$ I ran this to $t=4$ and got decent results, you may want to run longer. The above solution can be derived by letting $\partial_t\psi=0$ and solving the resulting ODE by normal techniques.

  • $\begingroup$ Is there anywhere that I might find a worked example of this, so I can write a test case for my code? $\endgroup$
    – jmansion
    Sep 30, 2014 at 5:41
  • $\begingroup$ @jmansion: I've updated my answer to include the missing subscript and an analytic solution for testing purposes. $\endgroup$
    – Kyle Kanos
    Sep 30, 2014 at 14:26
  • $\begingroup$ This has to be like the sixth or seventh time I've seen this bumped by community :( $\endgroup$
    – Kyle Kanos
    Dec 26, 2019 at 2:35
  • $\begingroup$ @KyleKanos Hopefully not again! $\endgroup$
    – tpg2114
    Dec 26, 2019 at 13:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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