PDA

View Full Version : von neumann stability analysis


deadringer
Apr28-07, 04:57 PM
We have a hyperbolic pde (in fact the 1d wave equation) with indep vars X, T
We use the central difference approximations for the second derivatives wrt X, T to get

[phi(Xn, Tj+1) -2phi(Xn, Tj) + phi(Xn, Tj+1)]/(dT^2) = [c^2][phi(Xn-1, Tj) -2phi(Xn, Tj) + phi(Xn=1, Tj)]/(dX^2)

where dX and dT are the changes in the independent variables between steps, c is the speed of the wave.

We are supposed to show that the scheme is stable only if cdT/dX = lambda<=1

We have to apply Von Neumann stability analysis. This usually means assuming a solution of the form phi = A(T)*exp(ikX)

Xn = ndX

therefore phi(Xn, Tj) = Aj*exp(ikndX)

The stability condition is that sigma = Aj+1/Aj and we want the modulus of sigma is less than or equal to one.

The first problem I have is that the above PDE finite difference equation has phi's with three different times in it; j-1,j,j+1. Therefore I tried using two differenct sigma factors sigma1 = Aj/Aj-1 and sigma2 = Aj+1/Aj
I do not think that you can assume that these are the same. I think that for stability you need the modulus of both sigmas is less than or equal to one.

the original finite difference equation gives us

sigma1 = 1/[2(lambda^2)cos(kdX) +2(1-(lambda^2)) - sigma2)

I set the modulus of this is less than or equal to one (call this stage A), giving two more equations:

-2(lambda^2)cos(kdX) - 2 + 2(lambda^2) + sigma2 <=1

this must be true for the largest values of the LHS (i.e when sigma2 = 1) , giving us Lambda^2 <= 1/(1-cos(kdX))

The RHS is largest when cos = 0, and also note that dX, dT, c are all positive by def'n, giving us the required condition. The only problem is we also get another equation from stage A:

2(lambda^2)cos(KdX) + 2 - 2(lambda^2) - sigma2 >=1

This must be true for the smallest values of the LHS i.e when sigma2 = 1
This implies that (lambda^2)(cos(kdX) - 1)>=0

We know that lambda^2>=0 which means that the other bracketed term must be greater than or equal to one. This is only satisfied for cos = 0
This extra condition was not mentioned in the question and doesn't look right. I'm not sure where I've gone wrong.

LeBrad
Apr29-07, 12:14 PM
From what I remember about this stuff, you plug the assumed solution \phi(x,t)=A(t)e^{ikx} into the difference equation, and then use the fact that x[n+m]-x[n]=m \Delta x to obtain an ordinary difference equation in t_j which can be solved exactly.

So in this case you would obtain
e^{ikx[n]}(A[j+1]-2A[j]+A[j-1]) = \frac{c^2 \Delta t^2}{\Delta x^2} A[j] (e^{ikx[n-1]} - 2e^{ikx[n]} + e^{ikx[n+1]})
which you can rearrange as

A[j+1] + (\alpha-2)A[j] + A[j-1] = 0


where
\alpha = \frac{c^2 \Delta t^2}{\Delta x^2}(2cos(k\Delta x)-2)
( or something similar to that).

Difference equations have solutions of the form A[j] = z^j, so plug that solution into the final difference equation above and you end up with a quadratic equation in z. Note that the solutions z depend on alpha. Then, since the original solution is now \phi = z^j e^{ikx}, you need z<=1 to prevent the solution from blowing up as j grows, so you need to pick alpha to give you z roots which are less than one.

So if you solve the quadratic equation for z as a function of alpha, you will get a constraint on alpha, and then using the definition of alpha, that should give you a constraint on c^2dt^2/dx^2 which is the one you're after.