[zdunol]

Hey guys,

I am working through Numerical Solution of Partial Differential Equations by K.W Morton and David Mayers. I'm doing forward and backward difference scheme for simple heat equation (u_t=u_xx) and authors say that in order to find solution with backward time difference (when we don't know u_n+1) we need Thomas algorithm. But to be honest I can see very little difference between backward and forward diff. scheme. Only thing is that backward scheme requires info about what is going on in the future, so if I say that for domain x=[0,1], t=[0,1] I got info about u=[:,1]=u_o I just go backward in time up to u[:,0] and everything is analogous to the forward scheme.
What am i missing guys?

[mprinkey]

The difference is huge. Forward Euler is only conditionally stable. In fact, you may only choose a timestep that is proportional to dx^2, For finely resolved meshes, this is crippling!

Backward Euler is unconditionally stable (for the heat equation anyway), but that Thomas algorithm is an answer for 1D problems. For 2D or 3D problems, there is no trivial O(N) algorithm to solve for future time. So, you get the option to choose timesteps to be whatever size you like, but you have to invest in the coding and solution time associated with a sparse linear equation solution.

From a mathematical level, you can look at forward Euler as a Talyor expansion around the current time. Backward Euler is a Taylor expansion about the future time. Both have the same formal truncation error, but by treating all future time values as unknowns, there is an implicit consistency to the results. Forward Euler has consistency too, but at the present time---nothing for the future time. So, if you choose too large of a timestep with forward Euler, you will get zigzagging nonsense and eventually NaN or INFs. If you do the same with backward Euler, you will just get the steady state solution (as time -> infinity) dictated by the boundary conditions.

[zdunol]

Thank you very much for the response, but i still do not understand the difference well. If I uderstand it corectly and forward method goes like this for heat equation
U(j,n+1)=U(j,n)+dt/dx^2*(U(j+1,n)-2*U(j,n)+U(j-1,n)) and backward method is like this:
U(j,n)=U(j,n+1)-dt/dx^2+(U(j+1,n+1)-2*U(j,n+1)+U(j-1,n+1)
then there is very little difference between the two methods, the differences in stability arises due to + vs - operation in there. I just don't know why Thomass algorithm is incorporated whereas in order to calculate stuff we just need info about n+1 step? ;c

[FMDenaro]

Originally Posted by zdunol

Thank you very much for the response, but i still do not understand the difference well. If I uderstand it corectly and forward method goes like this for heat equation
U(j,n+1)=U(j,n)+dt/dx^2*(U(j+1,n)-2*U(j,n)+U(j-1,n)) and backward method is like this:
U(j,n)=U(j,n+1)-dt/dx^2+(U(j+1,n+1)-2*U(j,n+1)+U(j-1,n+1)
then there is very little difference between the two methods, the differences in stability arises due to + vs - operation in there. I just don't know why Thomass algorithm is incorporated whereas in order to calculate stuff we just need info about n+1 step? ;c

Let me write the implicit scheme as:
U(j,n+1)-dt/dx^2+(U(j+1,n+1)-2*U(j,n+1)+U(j-1,n+1) = U(j,n)

That represents an algebric system in the unknown vector at time n+1:
A * U(n+1) = U(n)

the matrix A has only tri-diagonal entries and Thomas is well suited to solve it.

Note that, conversely, the explicit formulation leads to a decoupled system of equations, each one solved for a component U(j,n+1).

The two methods have the same accuracy when the local truncation error is analysed (first order in time and secondo order in space) but the stability constraints are totally different

[mprinkey]

You are missing the fundamental issue. You KNOW all of the U(j,n) values. You have those numbers in an array. For Foward Euler, the U[j,n+1] values are all EXPLICITLY computed from that array. That is one for/do loop and you have the value at the new time level. That is maybe 20 lines of code. Very easy. If you haven't yet, sit down and code it up. Fix the first and last element in the array as boundary conditions and don't update them. Set the rest of the array as some initial condition of your choice, and run it it with different timesteps. You will quickly see what we are talking about.

Now, look at Backward Euler. You don't have ANY of the U(j,n+1) values at all. You have U(j,n) values. And you have some algebraic relationships that look something like:

U(2,n+1) - 2*U(3,n+1) + U(4,n+1) = ....

for all of your values of j. Here j = 3 (the middle one). These equations couple the set of three points together. It is that set of algebraic equations (3 coefficients per equation and one equation for EACH internal value of j) that gets fed into the Thomas algorithm. That is a very different bit of code that you have to write. Not just one for/do loop creating an explicit update. You have to create the right hand side for the Thomas algorithm, and fill the coefficient matrix, and then solve it. I code probably do that in 100 lines of code, but it is much more complicated.

EDIT: You just beat me, Professor. 8)

[zdunol]

ok, so in implicit scheme i do not have the information about last step. That's what I suspected. Thank you all very much !