3D Diffusion Equation in MATLAB

[johnnyTransform]

Hi guys,

I have functioning MATLAB code for my solution of the 3D Diffusion equation (using a 3D Fourier transform and Crank-Nicolsen) that runs just from the command window and automatically plots the results. However, it seems like my solution just decays to zero regardless of what initial condition I use, or boundary conditions I implement. Could I toss this code to somebody so they can run it on their rig and let me know if they can diagnose it? I've been stuck on it for days and I'd much appreciate it.

If you'd like me to post it all in the thread using the code formatting tool, I can do that, but I'll take the suggestion first before completely diluting the post. Thanks guys

 

[johnnyTransform]

Just as a conjecture for those who are curious, I think it may have to do with the fact that I'm forwarding my Fourier coefficients in time using an explicit scheme. However, it's still an open question

EDIT: Sorry, I couldn't find the option to edit my original post. Not trying to bump this up
EDIT 2: Thank you mods, this subforum seems more appropriate for the question. Much appreciated!

 

[johnnyTransform]

Update: I think I may have figured out my universal decay problem. My solutions now display heat flux, and not just a static decay to zero. I'll explain what I'm thinking might be the trick for people who have similar problems in the future:

Starting with the 3D Diffsuion Eq.: $u_t=k(u_{xx}+u_{yy}+u_{zz})$, assume the solution has a Fourier expansion ($e^{-ipz},e^{-imy},e^{-inx}$).

Physical space: $u(x,y,z,t)$. Fourier space: $c(m,n,p,t)$.

So we have a continuous Fourier transform pair. Once you define your initial condition on a grid, that continuous Fourier transform pair becomes a discrete Fourier transform pair (since your solution becomes discrete once you define it on a grid). In the DFT pair, {n,m,p} is restricted from $(0,N-1)$, where N is the number of spatial points in your grid.

So, when you're forwarding your Fourier coefficients in time, and you need to loop through your spatial frequencies {n,m,p}, make sure to restrict your loops to positive values, so for example, in MATLAB:

for n=1:1:N
    for m=1:1:N
        for p=1:1:N

Not:

for n=-N/2:1:N/2
    for m=-N/2:1:N/2
        for p=-N/2:1:N/2

which would be appropriate if we were still dealing with our continuous Fourier transform pair, which allows for spatial frequencies {n,m,p} over $(-\infty,\infty)$. The continuous Fourier transform pair is what most people are familiar with, which explains why this would be a natural mistake to make. With this in mind, you can use fftn() and ifftn() to loop back and forth between Fourier space and physical space, while forwarding your Fourier coefficients in time using some type of numerical scheme (for example, Crank-Nicolsen). Hope this is helpful (and correct)!