1D diffusion equation via Jacobi (or other) method
Author
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,
Author
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?
Does this sound like a reasonable hypothesis?
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.
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.
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.
Author
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:
Implicit (Jacobi):
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:
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?
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.69814156157120500000You 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?
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.
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]
Author
Thanks for the response!
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'.
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?
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?
You are multiplying by the current time rather than the time step in your implicit method. I believe your code should be:
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.
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.
Author
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.
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.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
