A Numerically solving a transport equation

hunt_mat
Homework Helper
Messages
1,798
Reaction score
33
TL;DR Summary
I have a simple PDE that I want to solve:
[tex]\frac{\partial v}{\partial t}=-\frac{\partial v}{\partial x}[/tex] along with boundary conditions:
[tex]\frac{\partial v}{\partial x}\Bigg|_{x=0,1}=0[/tex]
I know this has a simple analytical solution, but that's not the point.
I'm using a ``downwind'' approximation for the spatial derivative:
\frac{\partial v}{\partial x}\approx -\frac{3}{2h}v_{j}+\frac{2}{h}v_{j-1}-\frac{1}{2h}v_{j-2}
I'm using the usual approximation for the time derivative, I get the following for a stencil:
v_{i+1,j}=\left(1+\frac{3\alpha}{2}\right)v_{i,j}-2\alpha v_{i,j-1}+\frac{\alpha}{2}v_{i,j-1}
where \alpha=\delta t/\delta x. To deal with the endpoints I evaluate the governing equation at the boundary points to get
\frac{\partial v}{\partial t}\Bigg|_{x=0,1}=0, and so this simply makes:
v_{i+1,1}=v(i,1),\quad v_{i+1,N}=v_{i,N}
To deal with the point at j=2, I write the BC at x=0 as
\frac{v_{i,1}-v_{i,0}}{\delta x}=0
This means that v_{i,0}=v_{i,1}. This can be used in the stencil at j=2 to yield
v_{i+1,2}=\left(1+\frac{3\alpha}{2}\right)v_{i,2}-\frac{3\alpha}{2}v_{i,1}
This completes the numerical model. It looks okay to me but it's unstable, and I don't see how. Can anyone suggest anything?
 
Physics news on Phys.org
What is the Courant Number you use? (##C = \frac{u\Delta t}{\Delta x}##) Since this is an explicit scheme, it should be lower than one. Also, you want to be 'upwinding' not 'downwinding', i.e. the point you calculate should be dependent on the points upwind from it, since that is what will be influencing that point.
 
My Courant number is 0.1. I should use:
\frac{\partial v}{\partial x}\Bigg|_{x=x_{j}}\approx-\frac{3}{2\delta x}v_{j}+\frac{2}{\delta x}v_{j+1}-\frac{1}{2\delta x}v_{j+1}
for my approximation?
 
Just as a side remark:
normally the counter ##i## would be used for x-direction, ##j## for the y-direction, ##k## for the z-direction and ##n## for time. Same for the ##\delta## sign. If you mean a small finite distance then use ##\Delta## instead.
It would be very convenient if you use the same counters since that makes everything you write down much more readable. I will definitely use it this way since otherwise I confuse myself.

I based my remark on your own remark in that you use the 'downwind' approximation. But it really depends on the sign of ##a## in ##\frac{\partial v}{\partial t} + a \frac{\partial v}{\partial x} = 0##. In your case ##a=1## and thus the upwind side is at the 'left' side of the query point (i.e. where x is lower). So if indeed ##x_{i-1} < x_{i}## then you are already using upwind, although you called it downwind.

Actually, it seems that you have a sign error in the equation below:
hunt_mat said:
I'm using a ``downwind'' approximation for the spatial derivative:
\frac{\partial v}{\partial x}\approx -\frac{3}{2h}v_{j}+\frac{2}{h}v_{j-1}-\frac{1}{2h}v_{j-2}

Since for a second order scheme it would look like this (which is the negative of what you got):
<br /> \frac{\partial v}{\partial x}\Bigg|_{x=x_{i}}\approx \frac{3v_{i} -4v_{i-1} +v_{i-2}}{2\Delta x}<br />Lastly, you now use a second order spatial derivative. If all fails, go back to a first order scheme since that is usually more robust.
 
Thread 'Direction Fields and Isoclines'
I sketched the isoclines for $$ m=-1,0,1,2 $$. Since both $$ \frac{dy}{dx} $$ and $$ D_{y} \frac{dy}{dx} $$ are continuous on the square region R defined by $$ -4\leq x \leq 4, -4 \leq y \leq 4 $$ the existence and uniqueness theorem guarantees that if we pick a point in the interior that lies on an isocline there will be a unique differentiable function (solution) passing through that point. I understand that a solution exists but I unsure how to actually sketch it. For example, consider a...
Back
Top