Problem with Crank-Nicholson for solving hyperbolic equations

[uqjzhang]

I am a noob so please let me know if here isn't the right place to post this.

Recently I am trying to solve hyperbolic equation m*dU/dt=k*d^2U/dx^2+q using Crank-Nicholson and finite element method. The final form of the solution is to compute the increment of the unknown at each time step through: (M/dt+0.5*K)dU=-KU(t)+q, U(t+dt)=U(t)+dU

The problem I am solving is quite simple: for a 1D bar x∈[0,1], initial U(x,t0)=U0=1, m=k=1. Boundary conditions are: U(1,t)=U0=1 (right end fixed); q(0.5,t)=-0.01 (point sink at center of the bar).

I read papers about time step control about this kind of problem. Obviously the time step can not be too small or it will cause oscillation around the sink location. But when I choose large time steps, I get this. (in this case, dt=3000s)
http://www.freeimagehosting.net/uploads/b541c099d4.jpg

I don't understand one thing here is: why the dU for these 4 time steps are so uneven? dU at odd time steps (e.g. 1,3) are much larger than even time steps (e.g. 2,4). If I increase the time step size, problem gets more severe.

If you have ideas about how this happens or how to eliminate this, I will be appreciated if you can share some light with me.

 

[AlephZero]

You can think of the finite difference scheme as dividing up the transfer of heat energy in the system in two ways.
1. Between two adjacent grid points along the bar.
2. Between two consecutive times, at the same grid point.

To get "good" results with the minimum amount of computer time, the two things need to be the same order of magnitude.

If the time steps are too small compared with the mesh size along the bar, you are wasting effort because the model of the temperature along the bar is too coarse to represent the details of how the solution changes with time.

if the time steps are too big, you get effects like your plots, because the correct temperature distribution is changing faster than your timesteps can represent it.

The practical "experimental" thing to do is just halve your time step, and keep halving it till the oscillations go away.

If you want to be more "scientific" about it, do some research on a non-dimensional quantity called the Fourier Number and how it relates to finite difference methods.

There is another quantity called the Biot number which is relevant for convection boundary conditions (but not for this problem).

 

[uqjzhang]

Thank you AlephZero. Your reply helps me a lot on understanding the problem itself. About Fourier number and Biot number, I heard them before and I definitely going to check them out.
I tried different theta values in Crank-Nicholson method, e.g. U(t+dt)=U(t)+theta*dU
It turns out that with bigger value of theta, I almost solve the problem! See the picture.
http://www.freeimagehosting.net/uploads/9ae065f8b7.jpg

The dU at the sink point of different time steps are almost even now. I am satisfied with the result, but I still do not understand how it works. It seems with bigger theta value, the solution gets more 'implicit' other than 'explicit', is that right?

 

[AlephZero]

By separating the variables, you can find general solutions of the PDE of the form
(A sin px + B cos px) e^(-qt)
where there is a equation connecting p and q, which involves matieral properties.
What that says is that a fixed point in space, the solution usually looks something like a decaying exponential function of time.

So back off a bit from the PDE and think about a first order ODE where the solution is an exponential decay:
x' + x = 0, x(0) = 1.
The solution is x = exp(-t).

You can solve this numerically with difference schemes like
x(h) = x(0) + h[(1-k)x'(0) + k x'(h)]
where k is a fixed constant.
If k = 0 this is the explict Euler forward difference scheme.
if 0 < k <= 1 it is an implicit scheme.
It is first-order accurate unless k = 1/2, when it is second-order accurate.

This behaves the same way as the time-integration part of Crank Nicholson but it is simpler to see what is going on.

When k = 1/2 we have
x(h) = x(0) + h([x'(0) + x'(h))/2
Using the ODE to eliminate x', and the boundary condition to eliminate x(0),
x(h) = 1 - h(1 +x(h))/2

x(h) = (1-h/2)/(1+h/2)

compared with the exact solution
x(h) =exp(-h).

If you plot graphs of these two functions, you see they are very similar for small values of h. But when h > 2, the numerical solution for x(h) becomes negative and as h tends to infintiy, the numerical solution tends to -1.

So if you take several timesteps for h > 2, the numerical solutions will be a "damped" oscillation that decays to zero, not the exact monotonic decreasing function.

If you do this analysis for k > 1/2, you will find the oscillations decay quicker, but the solution is less accurate for small values of h. For k = 1 there are no oscillations at all but the accuracy is poor unless you use very small time steps.

For h < 1/2, you will find that |x(h)| > 1 for large values of h. This means the method would give a solution that is a growing oscillation. In other words it is unstable unless h is sufficiently small, which is what "conditionaly stable" means, of course.

You can extend this approach to a general method for error and stability analysis of difference schemes.

 

[uqjzhang]

Thank you so much AlephZero, your comments helps a lot on these issues.