Hi,

I came across a problem in which I need to discretize "d4u/dx4" term on a NON-UNIFORM grid using 2nd order Central Difference. I tried to derive the form using general Taylor series (u[i+2], u[i+1], u[i], u[i-1], & u[i-2]) but I strucked in between. The difficulty is the following:

a*(u[i+2]-u[i]) + b*(u[i+1]-u[i]) + c*(u[i-1]-u[i]) + d*(u[i-2]-u[i]) = function(du/dx, d2u/dx2, d3u/dx3, .....)

On the right hand side after expansion,

coeff of term [d5u/dx5] = 0

coeff of term [d4u/dx4] = 1

coeff of term [d3u/dx3] = 0

coeff of term [d2u/dx2] = 0

Solve these 4 relations to get a,b,c, and d. However, we still have the term du/dx on the right hand side whose co-efficient is NON-ZERO. My questions is that what to do with this term?? Should I incorporate more points in the Taylor-series (e.g. u[i+3] etc.) or just discretize the "du/dx" using CD-2 and use it in the above expression to get d4u/dx4?

If any of you knows about this, please let me know. It will be a great help to me.

Thank you very much in anticipation,

-Chandra

Write down the discrete formula for q=d^2u/dx^2 on your nonuniform grid then apply the formula to calculate d^2q/dx^2. Elimate q from the discrete system and you have your result.

Hey Tom, Thank you very much for your reply. I will try to implement it.

-Chandra

OK, here's the stencil:

df_{i} = cl*f_{i+CL} + bl*f_{i+BL} + al*f_{i+AL} + A*f_{i} + ar*f_{i+AR} + br*f_{i+BR} + cr*f_{i+CR} +

If the stencil had been designed for a uniform grid then

CL = -3; BL = -2; AL = -1; AR = +1; BR = +2; CR = +3;

If you want d^4/dx^4 on a nonuniform grid, you generally lose one order over the symmetrical stencil. Here's the first order stencil:

A = 24/(AL AR BL BR)

bl = 24/ ((AL - BL) (AR - BL) BL (BL - BR))

al = 24/ (AL (AL - AR) (AL - BL) (AL - BR))

ar = -24/ ((AL - AR) AR (AR - BL) (AR - BR))

br = -24/ ((AL - BR) BR (-AR + BR) (-BL + BR))

Truncation Error terms

0 *(xi^0)

0 *(xi^1)

0 *(xi^2)

0 *(xi^3)

1 *(xi^4)

((AL + AR + BL + BR)/5) *(xi^5)

and here's the 3rd-order stencil:

A =

(24 (BL BR + BL CL + BR CL + BL CR + BR CR + CL CR + AR (BL + BR + CL + CR) + AL (AR + BL + BR + CL + CR)))

----------------------------------------------------------------------------------------------------------

(AL AR BL BR CL CR),

cl =

24 (BL BR + BL CR + BR CR + AR (BL + BR + CR) + AL (AR + BL + BR + CR)) -----------------------------------------------------------------------

(AL - CL) (AR - CL) (BL - CL) (BR - CL) CL (CL - CR)

bl =

24 (BR CL + BR CR + CL CR + AR (BR + CL + CR) + AL (AR + BR + CL + CR)) -----------------------------------------------------------------------

(AL - BL) (AR - BL) BL (BL - BR) (BL - CL) (BL - CR)

al =

24 (BR CL + BR CR + CL CR + BL (BR + CL + CR) + AR (BL + BR + CL + CR)) -----------------------------------------------------------------------

AL (AL - AR) (AL - BL) (AL - BR) (AL - CL) (AL - CR)

ar =

-24 (BR CL + BR CR + CL CR + BL (BR + CL + CR) + AL (BL + BR + CL + CR))

------------------------------------------------------------------------

(AL - AR) AR (AR - BL) (AR - BR) (AR - CL) (AR - CR)

br =

-24 (BL CL + BL CR + CL CR + AR (BL + CL + CR) + AL (AR + BL + CL + CR)) ------------------------------------------------------------------------

(AL - BR) BR (-AR + BR) (-BL + BR) (BR - CL) (BR - CR)

cr =

-24 (BL BR + BL CL + BR CL + AR (BL + BR + CL) + AL (AR + BL + BR + CL)) ------------------------------------------------------------------------

(AL - CR) CR (-AR + CR) (-BL + CR) (-BR + CR) (-CL + CR)

Truncation Error terms

0 *(xi^0)

0 *(xi^1)

0 *(xi^2)

0 *(xi^3)

1 *(xi^4)

0 *(xi^5)

0 *(xi^6)

(BL BR CL + BL BR CR + BL CL CR + BR CL CR + AR (CL CR + BR (CL + CR) + BL (BR + CL + CR)) + AL (BR CL + BR CR + CL CR + BL (BR + CL + CR) + AR (BL + BR + CL + CR)) ) / 210 *(xi^7)

So, there's your answer. It was done in Mathematica using Fourier analysis. Good luck! Don't ask me for a tutorial ...

Hey Runge_Kutta, thank you very much for your detailed reply. I am going through it. Also, if I need, I will be looking for some relevent book related to the details you mentioned...to avoid a tutorial from you

By the way, this is really going to be quite useful,

thank you,

-Chandra

Hello Runge_Kutta, well, i need one more favor from you. Could you please tell me the same thing for d2u/dx2 also. Also, if you have some link to these details, that would also be very useful. I tried to find out in the mathematica's documentation on internet but unfortunately I couldn't. I am doing DNS for turbulent flow and that is why I need all these.

Thank you very much in advance, your help will be too helpful to me indeed.

-Chandra

Third-order version of d^2/dx^2

===============================

2 (BL BR + AR (BL + BR) + AL (AR + BL + BR))

{{A -> --------------------------------------------,

AL AR BL BR

2 (AR BR + AL (AR + BR))

> bl -> --------------------------------,

(AL - BL) (AR - BL) BL (BL - BR)

2 (BL BR + AR (BL + BR))

> al -> --------------------------------,

AL (AL - AR) (AL - BL) (AL - BR)

2 (BL BR + AL (BL + BR))

> ar -> ---------------------------------,

(AL - AR) AR (-AR + BL) (AR - BR)

-2 (AR BL + AL (AR + BL))

> br -> ----------------------------------}}

(AL - BR) BR (-AR + BR) (-BL + BR)

Truncation Error terms

0 *(xi^0)

0 *(xi^1)

1 *(xi^2)

0 *(xi^3)

0 *(xi^4)

[AR BL BR + AL (BL BR + AR (BL + BR)) ]*(xi^5)/60

============================================

I've never seen a book with this stuff. I hope you can read the Mathematica dump after this software mangles it!

Hey, thank you very much indeed. This will be too useful for me.

-Chandra