Jump to content
  • Advertisement
Sign in to follow this  
dangerdaveCS

1D diffusion equation via Jacobi (or other) method

This topic is 3498 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

I've been struggling like hell trying to figure out how to implement the diffusion/heat equation. The explicit method is fine and easy, but I want to be able to use large time steps so I need to use an implicit method. Ultimately, I want to use the inhomogenous diffusion equation in 3D with something a little more advanced than the Jacobi method. For now however, I've gone back to the drawing board and am trying to implement the simplest 1D diffusion equation possible, using Jacobi iteration. That is, grid spacing = 1, diffusion coefficient = 1, and time step = 1. With a a simple 1D array of 10 elements, with the boundary fixed to 0, and the centre (5th) element initialised to 100, I've used an explicit 0.0001 time step scheme and integrated for (1/0.0001) iterations, i.e. to simulate 1 second of time. This gives a result of 30.851 at the centre (5th) element. No matter what I try I cant replicate that result using an iterative method. I've scoured google like mad, got books from the library, and I think I understand the general idea, but I get insane answers. All the material I've read go on and on about matrix this and that, and anything that deals with diffusion specifically *solves* the diffusion equation entirely, i.e. doesnt stop after 1 second of time. The closest I've got was by following GPU Gems 3 navier-stokes implementation and just taking the viscous diffusion part of that, but even that simply doesnt give much info or any examples - just states some unintelligable equation and goes on to show a few scattered nonsensical fragment programs. Anyway, I'm ranting. I would be extremely grateful if someone could start me off here, give me some pseudo code for this extremely simple example or soemthing, then I can build on it. Really tearing my hair out here. Cheers,

Share this post


Link to post
Share on other sites
Advertisement
Actually, I think I have a theory... It turns out that explicit and implicit schemes give similar values for small numbers of explicit-scheme time steps, but the solutions start to diverge from one another, ending up by 1 second simulated time being completely different. I'm *hoping* this is actually an inaccuracy in the explicit scheme due to many iterations of rounding errors or whatever. In which case I don't have a problem at all.

Does this sound like a reasonable hypothesis?

Share this post


Link to post
Share on other sites
Look into the Crank-Nicholson method. The Numerical Recipes books have a discussion of this method for solving diffusion equations.

Regarding large times in the standard explicit finite difference scheme: You can use "scale" instead of "time" to allow for a geometric progression of scales rather than an arithmetic progression of times. I have some information and implementations at my filters page.

Share this post


Link to post
Share on other sites
Dave's recommendation of the Crank-Nicholson method for the heat equation is good. It's unconditionally stable and gives results that are quite close to exact (for the 1D case, depending on grid density). If you want to understand in general how to choose a method that works for your problem, look into von Neumann stability analysis, which is a well understood and elegant method for characterizing the stability of different integration schemes (implicit and explicit) for linear PDE's.

Share this post


Link to post
Share on other sites
Thanks for the replies guys. Right now I'm just trying to get the most basic setup possible working. I think it's time for some code:

Explicit:

void TimeStep(){
for (int i = 1; i < GRID_WIDTH-1; i++)
fOneDSystemDelta = (fOneDSystem[i+1]-2*fOneDSystem+fOneDSystem[i-1]);
for (int i = 1; i < GRID_WIDTH-1; i++)
fOneDSystem += fOneDSystemDelta * dt;
}



Implicit (Jacobi):

void Iterative(){
double time = iNumTimeStepsGone * dt;
for (int j = 0; j < 1000; j++)
{
for (int i = 1; i < GRID_WIDTH-1; i++)
fOneDSystemDelta = fOneDSystemBuff;
for (int i = 1; i < GRID_WIDTH-1; i++)
fOneDSystemBuff = (fOneDSystem + time*(fOneDSystemDelta[i+1] + fOneDSystemDelta[i-1]))/(1+2*time);
}
for (int i = 1; i < GRID_WIDTH-1; i++)
fOneDSystem = fOneDSystemBuff;
}



And a comparison of the results, taken from the 5th element of the extremely simple 1D system mentioned in my first post. E marks Explicit, I marks implicit:


1: E99.79999999050051000000
1: I99.80059799753210800000
2: E99.60059998105801100000
2: I99.60238409242622000000
3: E99.40179797167222400000
3: I99.40534653294992000000
4: E99.20359196934286700000
4: I99.20947372946422100000
5: E99.00597998804445900000
5: I99.01475425156009200000
6: E98.80896004870120900000
6: I98.82117682525692700000
7: E98.61253017916200700000
7: I98.62873033026072300000
8: E98.41668841417549600000
8: I98.43740379728095000000
9: E98.22143279536524600000
9: I98.24718640540437800000
10: E98.02676137120499800000

...<snip>...

990: E31.02255186302058800000
990: I44.87886456555588000000
991: E31.00362195416289900000
991: I44.86069506535765300000
992: E30.98472424067113100000
992: I44.84254728179450200000
993: E30.96585864324026000000
993: I44.82442117089545500000
994: E30.94702508279727500000
994: I44.80631668881356000000
995: E30.92822348050043500000
995: I44.78823379182548300000
996: E30.90945375773852200000
996: I44.77017243633104500000
997: E30.89071583613010600000
997: I44.75213257885275200000
998: E30.87200963752280500000
998: I44.73411417603537600000
999: E30.85333508399254800000
999: I44.71611718464553300000
1000: E30.83469209784284400000
1000: I44.69814156157120500000



You can see the answers start similar but diverge, becoming very very different by 1000 time steps (= 1 second simulated time).

Is there something wrong with my implicit function above? Is the accuracy of one or the other method defective in some way?

Share this post


Link to post
Share on other sites
I believe your implementation for Iterative is not correct. First, inside your j-loop you need a statement "time += dt". I initialized "time = 0.0" outside the loop. With this change, I get nearly the same numbers you posted at j = 1000. So you probably have this correct in your own code.

Second, the calculation of fOneDSystemBuff does not look right to me. With an implicit method, you should have a linear system to solve. If U[n][j] is the diffused value at time step n and at spatial index j, then you should have the equation: U[n][j] = -time*U[n+1][j+1] + (1+2*time)*U[n+1][j] - time*U[n+1][j-1]. It appears as if you solved for U[n+1][j] = U[n][j] + time*(U[n+1][j+1] + U[n+1][j-1])/(1+2*time) and used this in the right-hand side of the assignment of fOneDSystemBuff. The problem is that you do not yet know U[n+1]
  • values since you are currently at time step n. The equation I mentioned represents a tridiagonal system of equations that you must solve.

    Share this post


    Link to post
    Share on other sites
    Thanks for the response!

    Quote:
    Original post by Dave Eberly
    I believe your implementation for Iterative is not correct. First, inside your j-loop you need a statement "time += dt". I initialized "time = 0.0" outside the loop. With this change, I get nearly the same numbers you posted at j = 1000. So you probably have this correct in your own code.


    What I'm doing with the Jacobi method is simply to compare the result with the evolution of the explicit method. dt is 0.001 and the explicit function is run for 1000 time steps to give 1 second simulated time. For comparison purposes I've run the implicit function to show what the result would be at each moment of simulated time, i.e. iNumTimeStepsGone*dt = 'amount of simulated time so far via explicit method'.


    Quote:
    Original post by Dave Eberly
    Second, the calculation of fOneDSystemBuff does not look right to me. With an implicit method, you should have a linear system to solve. If U[n][j] is the diffused value at time step n and at spatial index j, then you should have the equation: U[n][j] = -time*U[n+1][j+1] + (1+2*time)*U[n+1][j] - time*U[n+1][j-1]. It appears as if you solved for U[n+1][j] = U[n][j] + time*(U[n+1][j+1] + U[n+1][j-1])/(1+2*time) and used this in the right-hand side of the assignment of fOneDSystemBuff. The problem is that you do not yet know U[n+1]
  • values since you are currently at time step n. The equation I mentioned represents a tridiagonal system of equations that you must solve.


  • OK, this sounds like it could be zeroing in on the problem. I'm trying to use a convergent relaxation technique, i.e. jacobi method, so I can use U[n+1], it will just be a guess until it converges. Can you give me some code or link or explanation of how the Iterative() function *should* be written?


    Share this post


    Link to post
    Share on other sites
    You are multiplying by the current time rather than the time step in your implicit method. I believe your code should be:

    void Iterative(){
    double time = iNumTimeStepsGone * dt;
    for (int j = 0; j < 1000; j++)
    {
    for (int i = 1; i < GRID_WIDTH-1; i++)
    fOneDSystemDelta = fOneDSystemBuff;
    for (int i = 1; i < GRID_WIDTH-1; i++)
    fOneDSystemBuff = (fOneDSystem + dt*(fOneDSystemDelta[i+1] + fOneDSystemDelta[i-1]))/(1+2*dt);
    }
    for (int i = 1; i < GRID_WIDTH-1; i++)
    fOneDSystem = fOneDSystemBuff;
    }


    The two methods were diverging because they were giving solutions for different times.

    Also note that you are using Gauss-Seidel for your implicit method, not Jacobi. Using the Jacobi method to solve the tridiagonal system would make it identical to the explicit method. Both methods are guaranteed to converge since the matrix is diagonally dominant, but the Gauss-Seidel method should converge faster than the Jacobi method.

    Share this post


    Link to post
    Share on other sites
    Quote:
    Original post by Hexus
    You are multiplying by the current time rather than the time step in your implicit method. I believe your code should be:

    *** Source Snippet Removed ***

    The two methods were diverging because they were giving solutions for different times.

    Also note that you are using Gauss-Seidel for your implicit method, not Jacobi. Using the Jacobi method to solve the tridiagonal system would make it identical to the explicit method. Both methods are guaranteed to converge since the matrix is diagonally dominant, but the Gauss-Seidel method should converge faster than the Jacobi method.


    Thanks for the response!

    OK, my understanding of Gauss-Seidel is that it would be:


    void Iterative(){
    double time = iNumTimeStepsGone * dt;
    for (int j = 0; j < 1000; j++)
    {
    for (int i = 1; i < GRID_WIDTH-1; i++)
    fOneDSystemBuff = (fOneDSystem + time*(fOneDSystemBuff[i+1] + fOneDSystemBuff[i-1]))/(1+2*time);
    }
    for (int i = 1; i < GRID_WIDTH-1; i++)
    fOneDSystem = fOneDSystemBuff;
    }



    I.e. using the new value of fOneDSystemBuff immediately for subsequent grid points. Please, please do correct me if I'm wrong here.

    As for the time step, iNumTimeStepsGone is the number of time steps so far done via the explicit method. So if the explicit method was on time step 5, with dt=0.001, then there has been 0.005 seconds of time simulated via the explicit method. In order to check the value at this point against the implicit method I set the time step to iNumTimeStepsGone*dt (in this example =0.005), then use the implicit method. I've then put the results of this timestep-by-timestep comparison in my post above.

    I realize I've been a bit confusing with my code. In the end what I want to do is simulate 1 second of time, both via explicit and via implicit, and the results of implicit and explicit methods to match one-another.

    Share this post


    Link to post
    Share on other sites
    Quote:
    Original post by dangerdaveCS
    I'm trying to use a convergent relaxation technique, i.e. jacobi method, so I can use U[n+1], it will just be a guess until it converges.


    There are two iterations going on here. One is the iteration over time. The other iteration is associated with the "relaxation" in the Jacobi method to solve a linear system. I see only the time iteration in your code. The linear system you have to solve is tridiagonal, so you can do that in O(n) time (n is your grid width in this case). It is not clear to me why you would want the more general Jacobi method.

    Share this post


    Link to post
    Share on other sites
    Sign in to follow this  

    • Advertisement
    ×

    Important Information

    By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

    GameDev.net is your game development community. Create an account for your GameDev Portfolio and participate in the largest developer community in the games industry.

    Sign me up!