Discretising Elliptic PDE in cylindrical coordinates

[ognik]

Given an energy functional $ E=\int_{0}^{\infty} \,dr.r\left[\frac{1}{2}\left(\d{\phi}{r}\right)^2 - S.\phi\right] $
I am told that discretizing on a lattice r_i=ih (h=lattice size, i is i axis) leads to :

$ 2{r}_{i}{\phi}_{i} - {r}_{i+\frac{1}{2}}{\phi}_{i+1} - {r}_{i-\frac{1}{2}}{\phi}_{i-1} = {h}^{2}{r}_{i}{S}_{i} $

I am also told that that a lattice $ {r}_{i}=\left(i - \frac{1}{2}\right)h $ leads to the same result. And yes, I worked this out and it did lead to the above equation - but when I did the exercise for the lattice r_i=ih, I couldn't get the same result (see below):

From $ E=\int_{0}^{\infty} \,dr.r\left[\frac{1}{2}\left(\d{\phi}{r}\right)^2 - S.\phi\right] =\int_{0}^{\infty} \,dr\left[r\frac{1}{2}\left(\d{\phi}{r}\right)^2 - r.S.\phi\right] $

Discretising, $ E=\frac{1}{2h}\sum_{1}^{N} {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 - h \sum_{1}^{N} {r}_{i}{S}_{i}{\phi}_{i} $

Using the variation principle we ask that $ \pd{E}{{\phi}_{i}} = 0 $
$ So\: \pd{E}{{\phi}_{i}} = \pd{E}{{\phi}_{i}}[\frac{1}{2h}\sum_{1}^{N} {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 - h \sum_{1}^{N} {r}_{i}{S}_{i}{\phi}_{i}] = 0 $

$ \therefore \frac{1}{2h} \pd{E}{{\phi}_{i}}\sum_{1}^{N} {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 = h {r}_{i}{S}_{i} $
$ \therefore \frac{1}{2} \pd{E}{{\phi}_{i}} [{r}_{1}\left({\phi}_{1} - {\phi}_{0}\right)^2 + ... +\: {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 + \: {r}_{i+1}\left({\phi}_{i+1} - {\phi}_{i}\right)^2 +...+{r}_{N}\left({\phi}_{N} - {\phi}_{N-1}\right)^2 ]
= h^2 {r}_{i}{S}_{i} $

$ \therefore \left( {r}_{i} + {r}_{i+1}\right) {\phi}_{i} - {r}_{i}{\phi}_{i-1} - {r}_{i+1}{\phi}_{i+1} = h^2 {r}_{i}{S}_{i} $

This is close to the equation I am trying to get to, and I can use $ \frac{1}{2}\left({r}_{i} + {r}_{i+1}\right) = {r}_{i+\frac{1}{2}} $ to get
$ \therefore 2{r}_{i+\frac{1}{2}} {\phi}_{i} - {r}_{i}{\phi}_{i-1} - {r}_{i+1}{\phi}_{i+1} = h^2 {r}_{i}{S}_{i} $

Clearly if I subtract $\frac{1}{2}$ from all the indices I have got there for the left side (but not the right) - but (apart from the right side) how can I justify that? I would seem to be shifting the basis by $\frac{1}{2}$ - and with r from cylindrical coords, $ {r}_{i} \ne {r}_{i-\frac{1}{2}} $

 

[I like Serena]

Given an energy functional $ E=\int_{0}^{\infty} \,dr.r\left[\frac{1}{2}\left(\d{\phi}{r}\right)^2 - S.\phi\right] $
I am told that discretizing on a lattice r_i=ih (h=lattice size, i is i axis) leads to :

$ 2{r}_{i}{\phi}_{i} - {r}_{i+\frac{1}{2}}{\phi}_{i+1} - {r}_{i-\frac{1}{2}}{\phi}_{i-1} = {h}^{2}{r}_{i}{S}_{i} $

I am also told that that a lattice $ {r}_{i}=\left(i - \frac{1}{2}\right)h $ leads to the same result. And yes, I worked this out and it did lead to the above equation - but when I did the exercise for the lattice r_i=ih, I couldn't get the same result (see below):

From $ E=\int_{0}^{\infty} \,dr.r\left[\frac{1}{2}\left(\d{\phi}{r}\right)^2 - S.\phi\right] =\int_{0}^{\infty} \,dr\left[r\frac{1}{2}\left(\d{\phi}{r}\right)^2 - r.S.\phi\right] $

Discretising, $ E=\frac{1}{2h}\sum_{1}^{N} {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 - h \sum_{1}^{N} {r}_{i}{S}_{i}{\phi}_{i} $

Using the variation principle we ask that $ \pd{E}{{\phi}_{i}} = 0 $
$ So\: \pd{E}{{\phi}_{i}} = \pd{E}{{\phi}_{i}}[\frac{1}{2h}\sum_{1}^{N} {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 - h \sum_{1}^{N} {r}_{i}{S}_{i}{\phi}_{i}] = 0 $

$ \therefore \frac{1}{2h} \pd{E}{{\phi}_{i}}\sum_{1}^{N} {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 = h {r}_{i}{S}_{i} $
$ \therefore \frac{1}{2} \pd{E}{{\phi}_{i}} [{r}_{1}\left({\phi}_{1} - {\phi}_{0}\right)^2 + ... +\: {r}_{i}\left({\phi}_{i} - {\phi}_{i-1}\right)^2 + \: {r}_{i+1}\left({\phi}_{i+1} - {\phi}_{i}\right)^2 +...+{r}_{N}\left({\phi}_{N} - {\phi}_{N-1}\right)^2 ]
= h^2 {r}_{i}{S}_{i} $

$ \therefore \left( {r}_{i} + {r}_{i+1}\right) {\phi}_{i} - {r}_{i}{\phi}_{i-1} - {r}_{i+1}{\phi}_{i+1} = h^2 {r}_{i}{S}_{i} $

This is close to the equation I am trying to get to, and I can use $ \frac{1}{2}\left({r}_{i} + {r}_{i+1}\right) = {r}_{i+\frac{1}{2}} $ to get
$ \therefore 2{r}_{i+\frac{1}{2}} {\phi}_{i} - {r}_{i}{\phi}_{i-1} - {r}_{i+1}{\phi}_{i+1} = h^2 {r}_{i}{S}_{i} $

Clearly if I subtract $\frac{1}{2}$ from all the indices I have got there for the left side (but not the right) - but (apart from the right side) how can I justify that? I would seem to be shifting the basis by $\frac{1}{2}$ - and with r from cylindrical coords, $ {r}_{i} \ne {r}_{i-\frac{1}{2}} $

Hey Ognik,

What you write matches with the trapezoid rule.
You have evaluated a Right Rieman Sum instead of applying the Trapezoid Rule.
The difference matches more or less with a shift of a half grid size.

 

[ognik]

Hi I like Serena, thanks for the links.
Sorry, but I can't see what you mean, to me it looks like they use a 2-point difference for the differential and the integration becomes the summation approximation - could you please help me step through what you mean? Thanks.

 

[I like Serena]

Fair enough. It's not the trapezoid rule that makes the difference.
The problem is that the derivative in your discretization is half a grid size off.

The derivative at $r_i$ should be:
$$\pd \phi r(r_i) \approx \frac{\phi_{i+\frac 12} - \phi_{i-\frac 12}}{h} \approx \frac{\phi_{i+1} - \phi_{i-1}}{2h}$$

What if we discretize with:
$$ E=\frac{1}{8h}\sum_{1}^{N} {r}_{i}\left({\phi}_{i+1} - {\phi}_{i-1}\right)^2 - h \sum_{1}^{N} {r}_{i}{S}_{i}{\phi}_{i} $$

 

[ognik]

What happens I'm sorry to say, is new confusion :-) Why 1/8h please?
I went through the math for your suggestion, and got to:

$ \frac{1}{4h}\left[\left({r}_{i-1} + {r}_{i+1}\right){\phi}_{i} - {r}_{i-1}{\phi}_{i-2} - {r}_{i+1}{\phi}_{i+2}\right] $

I can use $ \frac{1}{2}\left({r}_{i-1} + {r}_{i+1}\right) = {r}_{i}$, but just can't see how to get to the required equation for the other 2 terms on the LHS from there.

I have checked back in the book and they do say to use 'the 2-point difference formula to approximate each first derivative at the points halfway between the lattice points'. They also say to use the trapezoidal rule for the integrals...

 

[ognik]

Hi anyone, would really appreciate some more help here, I know I must be missing something - probably basic - in my understanding of the way lattices work. I've been over this many times and just don't 'get it'. Thanks for all assistance :-)

 

[ognik]

Hi I like Serena (or anyone else who can help), where I am is that I understood and applied the below:
$$\pd \phi r(r_i) \approx \frac{\phi_{i+\frac 12} - \phi_{i-\frac 12}}{h} \approx \frac{\phi_{i+1} - \phi_{i-1}}{2h}$$

and ended up with: $ \frac{1}{2}\left[\left({r}_{i-1}+{r}_{i+1}\right){\phi}_{i} - {r}_{i+1}{\phi}_{i+2} - {r}_{i-1}{{\phi}_{i-2}}_{}\right]
={h}^{2}{r}_{i}{S}_{i} $

I used $ \frac{1}{4}\left({r}_{i-1}+{r}_{i+1}\right) ={r}_{i}$ for the first part, so far so good.
But now I have to argue that $ \frac{1}{2}\left({r}_{i+1}{\phi}_{i+2}+{r}_{i-1}{\phi}_{i-2}\right) ={r}_{i+\frac{1}{2}}{\phi}_{i+1}+{r}_{i-\frac{1}{2}}{\phi}_{i-1} $

Yes, $ \frac{1}{2}\left({r}_{i+1}+{r}_{i-1}\right) ={r}_{i+\frac{1}{2}}+{r}_{i-\frac{1}{2}} $ and $ \frac{1}{2}\left({\phi}_{i+2}+{\phi}_{i-2}\right) ={\phi}_{i+1}+{\phi}_{i-1} $

...which leaves me thinking that what I need to argue is an acceptable approximation, but I am unable to justify that mathematically?

(I tried sketching something like $ {\phi}_{i}\: against \: {r}_{i}$ , but I do not think that has any meaning (despite again appearing to agree with what I'm trying to prove))

 

[I like Serena]

and ended up with: $ \frac{1}{2}\left[\left({r}_{i-1}+{r}_{i+1}\right){\phi}_{i} - {r}_{i+1}{\phi}_{i+2} - {r}_{i-1}{{\phi}_{i-2}}_{}\right]
={h}^{2}{r}_{i}{S}_{i} $

Since the left hand side does not contain references to $h$, it is independent of the grid size.
So let's change the grid size of the left hand size.
Then we get:
$$ \frac{1}{2}\left[\left({r}_{i-\frac 12}+{r}_{i+\frac 12}\right){\phi}_{i} - {r}_{i+\frac 12}{\phi}_{i+1} - {r}_{i-\frac 12}{{\phi}_{i-1}}_{}\right]
={h}^{2}{r}_{i}{S}_{i}$$Btw, in my opinion, the expression you want to arrive at is flawed.
It refers to elements that are not in the grid.
What we really need, is an expression that only refers to elements that are in the grid.
The expression you ended up with, is the right one - it should not be forced into something it should not be.

 

[ognik]

Thanks ILS (although I confess to a little frustration with myself, I should be able to see things like the h argument for myself).

I am not surprised at your comments on the flaws, I have the impression that the book sometimes sacrifices accuracy to allow for a simpler illustration of the principles - in this case discretisation using the variation principle.

One question remains - we can change the grid size on the LHS, why can we at the same time leave the r_i on the RHS unchanged? My thought is, that is because it is a fixed point in the grid and so changing the grid size doesn't move the fixed point? (Probably not that well phrased :-))

(Finally, this topic covers the 2nd last chapter of this course, in the last chapter I have one more issue where I am stuck halfway through - subject is "Parabolic PDE algorithm" - if you have any spare time ...but completely understand if not)