Finite Differencing Error Help

[member 428835]

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 spatial 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!

 

[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 [itex]\frac{\Delta t}{\Delta x} < C[/itex]
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

 

[member 428835]

Thanks, I'll look into the von Neumann analysis.

so you increase your time step when your total run time tf increases?

if I make the spatial domain bigger, all else equal, my ##\Delta z## get's larger, so why would this be any different with time?

 

[Chestermiller]

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 = h²,
$$\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

 

[the_wolfman]

It looks like your trying to model an equation of the form:
[itex]\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}[/itex]

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

Where [itex] c [/itex] is a flow velocity and [itex]\alpha [/itex] 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 [itex]\delta t <\frac{2\alpha}{c^2} [/itex] and [itex]\delta t <\frac{\delta z^2}{2\alpha} [/itex].

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

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.

 

[member 428835]

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!

 

[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.

 

[Chestermiller]

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?

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

 

[member 428835]

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.

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.

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.
Chet

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 spatial 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 spatial nodes.

 

[Chestermiller]

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?

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

 

[member 428835]

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

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!

 

[Chestermiller]

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!

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

Chet

 

[member 428835]

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

Chet

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}.$$

Thanks again for your help!

 

[Chestermiller]

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}.$$

Thanks again for your help!

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

 

[member 428835]

Are you experiencing instability with the values of Δt and Δx that you are currently using?
Chet

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 spatial 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?

 

[Chestermiller]

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 spatial nodes.

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

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).

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?

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

Chet

 

[member 428835]

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)##.

 

[Chestermiller]

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)##.

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

Chet

 

[member 428835]

Yes, this is correct. What are you thinking about this so far?

 

[Chestermiller]

Yes, this is correct. What are you thinking about this so far?

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

 

[member 428835]

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}$$
Chet

I can show this, although I get
$$\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}$$

Basically we have $$\frac{\partial y}{\partial t} = \frac{\partial y}{\partial t}+\frac{\partial y}{\partial ZL}\frac{\partial ZL}{\partial t} = \frac{\partial y}{\partial t} + \frac{Z}{L}\frac{\partial y}{\partial Z}\frac{\partial L}{\partial t}$$ where the ##Z/L## can come out of the derivatives since ##y## does not depend on ##L## and ##L## does not depend on ##Z##. This is the multivariable chain rule. Since ##y## does not depend on ##L## we have ##y^{3/2}_{zz} = y^{3/2}_{ZZ}/L^2##.

The only discrepancy I have is with the negative sign. How do you recommend we proceed?

 

[Chestermiller]

$$dy=\left(\frac{\partial y}{\partial t}\right)_Zdt+\left(\frac{\partial y}{\partial Z}\right)_tdZ$$
$$\left(\frac{\partial y}{\partial t}\right)_z=\left(\frac{\partial y}{\partial t}\right)_Z+\left(\frac{\partial y}{\partial Z}\right)_t \left(\frac{\partial Z}{\partial t}\right)_z$$
$$\left(\frac{\partial Z}{\partial t}\right)_z=-\frac{z}{L^2}\frac{dL}{dt}=-\frac{Z}{L}\frac{dL}{dt}$$
So,
$$\left(\frac{\partial y}{\partial t}\right)_z=\left(\frac{\partial y}{\partial t}\right)_Z-Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}$$
So,
$$\left(\frac{\partial y}{\partial t}\right)_z=\left(\frac{\partial y}{\partial t}\right)_Z-Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}=\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$
So,
$$\left(\frac{\partial y}{\partial t}\right)_Z=Z\left(\frac{\partial y}{\partial Z}\right)_t\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2(y^{3/2})}{\partial Z^2}$$

 

[member 428835]

Cool, I'm on board here but what are the subscripts for?

 

[Chestermiller]

Cool, I'm on board here but what are the subscripts for?

The subscripts are standard notation when doing partial differentiation for indicating what parameters are being held constant.

Now, before discussing how to discretize the partial differential equation, I want to immediately derive an equation for determining how L is varying as a function of time. We accomplish this by working with the differential equation at Z = 1. Right at Z = 1, y = 0. If we substitute this into the differential equation, we obtain:
$$\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}=0\tag{1}$$
If we solve this equation for dL/dt, we obtain 0/0, so we need to evaluate the limits of the spatial terms very carefully to find out what the correct limit of the ratio is. Since both y and its derivative dy/dZ are zero at Z = 1, y must be varying as follows in close proximity to the Z = 1:
$$y=k(1-Z)^2$$
where, from the Taylor series expansion,
$$k=\frac{1}{2}\left(\frac{\partial^2 y}{\partial Z^2}\right)_{Z=1}$$
Therefore, in close proximity of Z = 1,
$$y^{3/2}=k^{3/2}(1-Z)^3$$
$$\frac{\partial y}{\partial Z}=-2k(1-Z)$$
$$\frac{\partial ^2(y^{3/2})}{\partial Z^2}=6k^{3/2}(1-Z)$$
Substituting these relationships into Eqn. 1 yields:
$$-2k(1-Z)\frac{1}{L}\frac{dL}{dt}+6k^{3/2}(1-Z)\frac{2}{3L^2}=0$$
Note that the (1-Z)'s can be cancelled, and we can now solve for dL/dt:
$$L\frac{dL}{dt}=2\sqrt{k}$$
and
$$\frac{dL^2}{dt}=4\sqrt{k}=2\sqrt{2\left(\frac{\partial^2 y}{\partial Z^2}\right)_{Z=1}}$$
This equation is subject to the initial condition L = L₀ @ t = 0.
These equations give us what we need to establish ##\frac{1}{L}\frac{dL}{dt}## in the differential equation, based on the second partial of y with respect to Z at Z = 1.

Chet

 

[member 428835]

Wow, that's good stuff Chet. I have a few questions though.

Since both y and its derivative dy/dZ are zero at Z = 1

How do you know dy/dZ is zero at Z = 1? I understand y = 0 at Z = 1.

where, from the Taylor series expansion,
$$k=\frac{1}{2}\left(\frac{\partial^2 y}{\partial Z^2}\right)_{Z=1}$$

How are you arriving at the above taylor series? Also you mentioned that the final solution is predicated on knowing this second derivative, which is possibly a function of time. How do you recommend we find this?

$$\frac{dL^2}{dt}=4\sqrt{k}=2\sqrt{2\left(\frac{\partial^2 y}{\partial Z^2}\right)_{Z=1}}$$
This equation is subject to the initial condition L = L₀ @ t = 0.

I also am missing why we need the above information regarding ##L^2##. Can you help?

Thanks so much for your patience!

 

[Chestermiller]

Wow, that's good stuff Chet. I have a few questions though.How do you know dy/dZ is zero at Z = 1? I understand y = 0 at Z = 1.

Your own boundary condition is dy/dz = 0 at z = L.

How are you arriving at the above taylor series?

$$y (Z) = y(1) + y'(1)(Z-1) + \frac{y''(1)}{2}(Z-1)^2=y(1) - y'(1)(1-Z) + \frac{y''(1)}{2}(1-Z)^2$$

Also you mentioned that the final solution is predicated on knowing this second derivative, which is possibly a function of time. How do you recommend we find this?

No problem. We'll get to that in a little while. We need to use a second order finite difference approximation for y'' at Z = 1 and time t=mΔt.

I also am missing why we need the above information regarding ##L^2##. Can you help?

We need to know L(t) by integrating with respect to time to get L²(t).

 

[member 428835]

Your own boundary condition is dy/dz = 0 at z = L.

In post 17 I have ##h_z(0,t) = 0##, which I believe implies ##h_Z(0,t)=0##, not necessarily at ##Z=1##. If this is true, the taylor expansion is complicated. What do you think?

Other than this I'm with you, thanks for the great explanations!

 

[Chestermiller]

In post 17 I have ##h_z(0,t) = 0##, which I believe implies ##h_Z(0,t)=0##, not necessarily at ##Z=1##. If this is true, the taylor expansion is complicated. What do you think?

In line with our transformation, at any time t,
$$\frac{\partial y}{\partial z}=\frac{1}{L(t)}\frac{\partial y}{\partial Z}$$
This equation applies at all locations of z (and corresponding values of Z=z/L, including z = L, Z = 1).
Chet

 

[Chestermiller]

In post 17 I have ##h_z(0,t) = 0##, which I believe implies ##h_Z(0,t)=0##, not necessarily at ##Z=1##. If this is true, the taylor expansion is complicated. What do you think?

Other than this I'm with you, thanks for the great explanations!

Oh. I had always thought that you had specified dy/dz = 0 at z = L. Checking back, I now see that you didn't. However, in order to conserve volume, we can show that dy/dz does have to be zero at z = L.

Chet

 

[Chestermiller]

Follow up to posts #18 and 19: Are you sure that the volume is the integral of h² (and not h)? The reason I'm asking is that this was what led to my conclusion that y' = 0 at z=L. If y' is not equal to zero at z = L, getting the time derivative of L is much simpler, and the limiting process I outlined in the last few posts is neither needed nor warranted.

Chet

 

[member 428835]

Follow up to posts #18 and 19: Are you sure that the volume is the integral of h² (and not h)?
Chet

Yes, I am sure that the volume can be represented as ##F \int_{-L}^L h^2 dz## where ##F## is a known constant. How were you thinking to proceed?

 

[Chestermiller]

Yes, I am sure that the volume can be represented as ##F \int_{-L}^L h^2 dz## where ##F## is a known constant. How were you thinking to proceed?

I have something definitely in mind. Let me get you started. What do you get if you take the derivative of V with respect to time, applying the Leibnitz rule on the right hand side of the equation?

Chet

 

[member 428835]

What do you get if you take the derivative of V with respect to time, applying the Leibnitz rule on the right hand side of the equation?
Chet

I assume ##V## is volume? Then we have $$\frac{\partial V}{\partial t} = \frac{\partial}{\partial t} \int_{-L}^L h^2 dz = \int_{-L}^L 2h h_t \, dz + 2h^2(L,t) \frac{dL}{dt}$$
since ##h## is even with respect to ##z##. Do you agree?

 

[Chestermiller]

I assume ##V## is volume? Then we have $$\frac{\partial V}{\partial t} = \frac{\partial}{\partial t} \int_{-L}^L h^2 dz = \int_{-L}^L 2h h_t \, dz + 2h^2(L,t) \frac{dL}{dt}$$
since ##h## is even with respect to ##z##. Do you agree?

Yes. Now h(L,t) = 0, so
$$\frac{dV}{dt}= \frac{\partial}{\partial t} \left(\int_{-L}^L h^2 dz\right) = \int_{-L}^L\left(\frac{\partial h^2}{\partial t}\right) dz $$
But, if we combine this with the differential equation, we obtain:
$$\frac{dV}{dt}= \int_{-L}^L\left(\frac{\partial h^2}{\partial t}\right) dz=\frac{4}{3}\left(\frac{\partial h^3}{\partial z}\right)_{z=L}=0 $$
Do you agree?

Chet

 

[member 428835]

Do you agree?

I do, since ##h## is even we evaluate from ##[0,L]## and multiply by ##2##. Then we are left with the expression you wrote since ##h_z(0,t)=0##. The last equality holds since volume is conserved. This is getting fun; what's the next plan? Obviously this qualifies your taylor expansion and post 24.