numerical stability due to cross derivatives
 

Dear friends:

I got a problem in my FDM method for surface wave model in 2D horizontal domain that I confirm instabilities are due to the treatment of cross derivatives. If any of you've seen any analysis of similar case, please tell me. More specifically I have a system like:

[u + u_xx + u_xy + u_yy + v_xx +v_xy + v_yy]_t = F(x,y,t)

[v + u_xx + u_xy + u_yy + v_xx +v_xy + v_yy]_t = G(x,y,t)

the way I solve it is like this:

Let U = u+u_xx, V= v+v_yy, then the equation becomes:

U_t = F(x,y,t) - u_xyt-u_yyt - v_xxt - v_xyt- v_yyt = F1(x,y,t)

and

V_t = G(x,y,t)-v_xyt - v_xxt - u_xxt - u_xyt - u_yyt = G1(x,y,t)

Then I use Adams projection in time to predict U and V for next time step (n+1) based on values obtained at n-2, n-1, n steps. And solve U = u+u_xx for u at n+1, also solve V= v+v_yy for v at n+1 step.

Problem being whenever I plug in the cross derivative terms such as u_xyt and v_xyt in the RHS of the above equation, then the system becomes unstable. If you have seen any analysis about this kind of problem, please let me know, many thanks,

Wen

Re: numerical stability due to cross derivatives
 

I have never seen anything like this. I am just speculating here, but following your own ideas, why not try solving the following system

(1) U = u + u_xx + u_xy + u_yy + v_xx +v_xy + v_yy

(2) V = v + u_xx + u_xy + u_yy + v_xx +v_xy + v_yy

(3) U_t = F(x,y,t)

(4) V_t = G(x,y,t)

First update U, V from t to t + Î" t using (3), (4). Then solve for u and v using (1), (2).

Re: numerical stability due to cross derivatives
 

Praveen:

the u_xx + u_xy + u_yy + v_xx +v_xy + v_yy part in the RHS of my equation are called dispersive terms for intermediate depth water waves. They goven the dispersion of waves (long waves go fast, short waves go slower such that the waves are separated over a distance of proporgation, the term "dispersive" is borrowed from light through a triganular shape glass which will disperse the light components to be a rainbow pattern)

I have tried moving everything to the right except u and v, the thing is that the F term and G term include convections such as u*u_x +v*u_y in F and u*v_x + v*v_y in G. These terms make the waves shoal and pitch forward and break, (shocks appear). Without solving the dispersion system U=u+u_xx , which is very simple dissipation type of equation, the solution will be easily unstable when shock forms. The occurance of u_xx basically smooth the shock to a handle-able degree.

You might be right, that I may should try it. The reason I didn't try is that in the nonlinear term F, and G, I did use an iteration and kind of don't want to solve the resulting equation in your suggested (1)(2) using iteration again. Now I'm solving U=u+u_xx and V= v+v_yy by tri-diagonal solver. If I solve U=u_xx+u_xy+u_yy ..., i would have to use iteration for a complex domain with curving coastlines and islands etc. Anyway, thanks for the suggestion, and that's probably what I have to do and sacrifice more computation time and should be much better with a stable code than an unstable code that still sort of rely on filtering.

My model is pretty available, but still needs a lot of work, you may want to take a look at http://www.coastal.udel.edu/~longmtm/tmp

I will report when I finally solve this problem.

Wen