Finite Differencing Error Help

1. Sep 18, 2015

joshmccraney

hi PF! I've attached pictures to help you all see what is happening. Basically, I am running a forward time-centered space finite difference scheme, which is $${h_i^{j+1} \approx \left[ h_i^j \left( \frac{h_{i+1}^j-2h_i^j+h_{i-1}^j}{\Delta z^2} \right) + \frac{1}{2}\left( \frac{h_{i+1}^{j}-h_{i-1}^j}{ \Delta z} \right)^2 \right]\Delta t+h_i^j}.$$ The system runs fine as long as my final run time $tf$ is low. Since $\Delta t$ is defined as $tf/M^2$ where $M$ is the number of spacial nodes (determined by scaling, which I know is accurate).

At any rate, the following pictures are from the same number of time iterations, but different the divergent picture has a much larger $tf$ value. Any idea why one diverges and the other maintains equilibrium? I know $tf$ is the reason, but why?

Thanks so much!

2. Sep 18, 2015

bigfooted

so you increase your time step when your total run time tf increases? That doesn't sound right, as most numerical methods for such PDE's have a CFL stability criterion of the form $\frac{\Delta t}{\Delta x} < C$
You can do a von Neumann analysis for your convection-diffusion problem and check the necessary and sufficient condition for stability:
https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis

3. Sep 18, 2015

joshmccraney

Thanks, I'll look into the von Neumann analysis.
if I make the spacial domain bigger, all else equal, my $\Delta z$ get's larger, so why would this be any different with time?

4. Sep 18, 2015

Staff: Mentor

As best I can tell, the original differential equation was:
$$\frac{dh}{dt}=h\frac{d^2h}{dz^2}+2\left(\frac{dh}{dz}\right)^2$$
Is there a reason you didn't do the following:
$$h\frac{dh}{dt}=h^2\frac{d^2h}{dz^2}+2h\left(\frac{dh}{dz}\right)^2=\frac{d}{dz}\left(h^2\frac{dh}{dz}\right)$$
So,
$$\frac{d(h^2)}{dt}=\frac{2}{3}\frac{d^2(h^3)}{dz^2}$$
So, if we define y = h2,
$$\frac{dy}{dt}=\frac{2}{3}\frac{d^2}{dz^2}y^{3/2}$$

Wouldn't this be easier to integrate?

Ooop. I forgot to use partials, but you get the idea.

Chet

5. Sep 18, 2015

the_wolfman

It looks like your trying to model an equation of the form:
$\frac{\partial h}{\partial t}-\frac{1}{2}\frac{\partial h}{\partial z}\frac{\partial h}{\partial z}=+h\frac{\partial^2 h}{\partial z^2}$

This looks like an ugly nonlinear form of the advection-diffusion equation:
$\frac{\partial h}{\partial t}+c\frac{\partial h}{\partial z}=\alpha\frac{\partial^2 h}{\partial z^2}$

Where $c$ is a flow velocity and $\alpha$ is a diffusivity of some sort. Using a explicit finite difference scheme to advance the linear advection diffusion equation, Von Neuman analysis shows that numerical stability requires $\delta t <\frac{2\alpha}{c^2}$ and $\delta t <\frac{\delta z^2}{2\alpha}$.

The analysis of the nonlinear problem is tricky, but we can come up with crude conditions for stability by recognizing that $c \sim -\frac{1}{2}\frac{\partial h}{\partial z}$ and $\alpha \sim h$. Using these with the above stability condition we get that
$\delta t < \frac{\delta z^2}{2h_{max}}$.

Again this is only a crude analysis. If you use a time step larger than the one I suggest, then it's pretty likely that you will get a numerical instability. However, using a time-step smaller than what is derived does not guarantee stability.

6. Sep 18, 2015

joshmccraney

Thank you both for replying, and I see Chet changed his signature favoring the Wolverines.

the_wolfman, I'm just learning about finite difference techniques, but I think von Neumann is only valid for linear problems since it implies a solution has separated variables. I have read, however, about linearizing the PDE and then applying von Neumann analysis. Does this sound right to you, or anyone?

And Chet, of course the PDE you proposed is the correct one I'm dealing with, and your calculus trick give an impressive result! I didn't re-write it because I simply didn't see it. So if I were to discretize the problem you have posed would I get something like this

$$\frac{y_i^{j+1}-y_i^j}{\Delta t} = \frac{2}{3}\frac{(y_{i+1}^j)^{3/2} - 2 (y_i^j)^{3/2} + (y^j_{i-1})^{3/2}}{\Delta y ^2}.$$

But this elegant rewritten form remains non-linear. Would you advise me, when checking for stability, to linearize it (which I've never done but am reading about it) and the use von Neumann or do you have something else up your sleeve?

Thank you both so much!

7. Sep 18, 2015

the_wolfman

It is Friday afternoon and my brain is not fully functional, so forgive me if I'm wrong. But I think Chet's approximation of your PDE is off by a coefficent. He has a factor 2 when I think it should be a half. I bring this up because this factor matters in converting the RHS of the equation into a total derrivative. Either way it's worth double checking.

Your right in that you often cannot apply von Neumann to nonlinear equations. One way to get a stability condition is to linearize the equations. However this will not guarantee stability. It we give you a rule of thumb. In using the result from the advection diffusion equation I kind of did this, but I miss the effect of two terms(?). My result does show why the stability should depend on both the time step and grid spacing.

If you do a more formal linearization I'd be curious as to what stability conditions you get.

Last edited: Sep 18, 2015
8. Sep 18, 2015

Staff: Mentor

First, why don't you try it to see if it works any better than the original scheme. I can think of lots of ways of improving the stability if you are not then happy.

Of course, the obvious stability parameter is $\frac{Δt}{(Δx)^2}$, but that's also going to depend on the entire profile at a given time.

Chet

9. Sep 18, 2015

joshmccraney

Chet's is correct. If you take the 1/2 in the equation i posted and change it to $2/2^2$ and then bring the $2^2$ into the fraction you get a $2 \Delta z$ in the denominator, which is what you want since this is a central difference scheme. Thanks for looking into it! Also, I'll post what stability conditions I arrive at.

I ran the y scheme you presented, and the results were bad. The issue is I think in the initial condition, which is in green above. Since values are less then 1 and lifted to a power greater than 1, they are all being dwarfed very small. As a result of just a few time loops I'm getting y values pretty much all zeros. Any ideas why?

But I am getting the original scheme I posted to run with larger values of $tf$ by heuristically changing the number of time loops by factor's of 10 without changing the number of spacial nodes. I'm mainly curious as to why increasing the time loops allows higher values of $tf$ to work. Any ideas or direction? Again, $dt :=tf/M^2$ where $M^2$ is the number of time iterations and $M$ is the number of spacial nodes.

10. Sep 18, 2015

Staff: Mentor

Try the calculation with an exponent of 1 on the y's just to make sure that the scheme is being implemented correctly. Then it's just a straightforward linear parabolic problem (diffusion equation).

Chet

11. Sep 18, 2015

joshmccraney

Dang, that's a smart approach (and intuitive too)! Sadly I have guests coming over tonight so I have stop working, but I'll pick back up tomorrow. Yes, there was a typo on my end which resulted in the code flatlining. It is doing the right thing now, so I'll transform back tomorrow and let you know my results.

On a side note, the PDE you wrote above (which is correct) models capillary flow, $h$ being the height and $z$ being downstream location, thus $h(ztip,t)=0$ is a boundary condition. It can be verified that at the tip, $z=ztip$, a linear solution in $z$ and $t$ solves the system, that is $h = a z + b t + c$ for constants $a$, $b$, and $c$. The entire code that I am using relies on stretching the domain on each end by using the previous two points and making a linear extrapolation. Would I simply make a quadratic extrapolation using the last two points since $y = h^2 = c^2 + 2 a c t + a^2 t^2 + 2 b c x + 2 a b t x + b^2 x^2$ is a quadratic?

If this is confusing let me know and I will make changes. Thanks so much!

12. Sep 19, 2015

Staff: Mentor

Yes, it's confusing, at least to me. I have no idea what it is saying of what the physical situation is.

Chet

13. Sep 21, 2015

joshmccraney

Sorry, I guess all I'm trying to say is the domain actually stretches. At any rate, how would you go about checking stability for

$$\frac{y_i^{j+1}-y_i^j}{\Delta t} = \frac{2}{3}\frac{(y_{i+1}^j)^{3/2} - 2 (y_i^j)^{3/2} + (y^j_{i-1})^{3/2}}{\Delta y ^2}.$$

14. Sep 21, 2015

Staff: Mentor

I would run the calculation over and over again with increasing values of Δt until the calculation went unstable. Are you experiencing instability with the values of Δt and Δx that you are currently using?

Chet

15. Sep 21, 2015

joshmccraney

No, the code works with the same values, but I've still managed to crash it. This occurs when $tf$ is sufficiently high and the number of time nodes is the square of spacial nodes.

Also, what is the appropriate way to perform finite difference scheme on the tip point? Would you create a ghost node or use a backward space formula?

16. Sep 21, 2015

Staff: Mentor

Better be careful. This could be equivalent to increasing Δt/(Δx)2.

Incidentally, I've worked out a scheme that will be much more stable, but I am reluctant to show it if your existing scheme gives you adequate performance (i.e., gives you stable solutions to the cases you are interested in).

I don't follow. Please state the differential equation and boundary conditions.

Chet

17. Sep 21, 2015

joshmccraney

The differential equation is $$h_t = 2(h_z)^2+h h_{zz}$$ The boundary condition is $h(L(t),t) = 0$ and $h_z(0,t)=0$. Additionally, the volume (which is a function of the area under the finite difference curve) is constant. What I have been doing is applying a symmetric initial condition like the quadratic i posted above. Then to find the location of $L(t)$ I use the previous two points to make a linear line. Where the height of that line equal zero is the new $L(t)$.

18. Sep 21, 2015

Staff: Mentor

So the constant volume is equal to:
$$V=\int_{-L(t)}^{+L(t)}h^2(t,z)dz$$
Is that correct?

Chet

19. Sep 21, 2015

joshmccraney

20. Sep 21, 2015

Staff: Mentor

I have most of the details worked out, but it might be a little beyond your experience, so I'm going to start slowly.

The first step is to transform the coordinates from t and z to t and Z = z / L(t). See if you show that, when you make this transformation, the differential equation transforms into:
$$\frac{\partial y}{\partial t}=Z\left(\frac{\partial y}{\partial Z}\right)\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$
The boundary conditions are now applied at Z = -1 and Z = +1.

Notice that the boundary conditions are very nice now, but now the complexity is transferred to the differential equation. And now, in order to apply the differential equation, you somehow need have an equation for dL/dt.

Chet