Why Does Increasing Final Run Time Cause Divergence in Finite Differencing?

[member 428835]

Ok, I think I corrected all the errors you pointed out. How does this attached version look?

As an aside, do we really need (12) through (16)? They seem to not really be used.

 

[Chestermiller]

You're awesome!

Not really. Just lots of experience.

 

[Chestermiller]

Ok, I think I corrected all the errors you pointed out. How does this attached version look?

As an aside, do we really need (12) through (16)? They seem to not really be used.

Yes. I would lose 12 - 16.

There's a typo in Eqn. 29: the y should be under the radical in the numerator.

Please show the derivation of Eqns. 30 and 34-36. I'm not able to follow in detail. I just want to make sure that it was done correctly. For example, I think there should be a 1/2 in Eqn. 30 instead of the factor of 2; and, in Eqn. 34, I feel like there should be a 4 in the denominator. Also, in Eqn. 36, there should be a ζ in the denominator rather than a ξ.

Chet

 

[member 428835]

How does this look?

(23) is evaluated at $t$ and yet requires $y$ to be evaluated at $t$ as well. How would you deal with this?

I also got rid of zeta and xi.

 

[Chestermiller]

How does this look?

I also got rid of zeta and xi.

I found a couple of errors in the new version:
Eqns. 17,19, 22,23,33,34,and 35 should have a 4 instead of a 2.

In Eqn. 32, you are missing an additional factor of 2 in the denominator of the first term in brackets. This is because ∂y/∂Z is evaluated numerically using a central difference over 2 grid intervals.

(23) is evaluated at $t$ and yet requires $y$ to be evaluated at $t$ as well. How would you deal with this?

In Eqns. 33-34, you need to replace the j's with (j+1) and you need to replace the (j-1)'s with j. And in Eqn. 35, you need to replace the (j-1) with j (so you have j on both sides). As far as the L_j² in the denominator of Eqn. 32 is concerned, you have already calculated that from the previous time step (Or, to the same degree of accuracy, you could also use L_j+1² or the average of the two).

 

[member 428835]

Thanks for pointing these errors out. I agree with all of your corrections, and believe this is now error free. Thanks for all the help!

How would I go about an error analysis though? It's non-linear so I don't think a von Neumann approach is valid since separation of variables doesn't work.

I'll start coding this too (in MatLab) and post a few pictures of what I'm getting in the .pdf for you're consideration.

 

[Chestermiller]

Thanks for pointing these errors out. I agree with all of your corrections, and believe this is now error free. Thanks for all the help!

How would I go about an error analysis though? It's non-linear so I don't think a von Neumann approach is valid since separation of variables doesn't work.

I'll start coding this too (in MatLab) and post a few pictures of what I'm getting in the .pdf for you're consideration.

Are you talking about error analysis, or are you talking about stability?

Chet

 

[member 428835]

Are you talking about error analysis, or are you talking about stability?

Chet

Sorry, I'm talking about stability.

 

[Chestermiller]

Sorry, I'm talking about stability.

I don't know of any methods of analyzing the stability of time-dependent schemes in non-linear problems, but, experience has shown that stability is guaranteed if you go to fully implicit (e.g., Backward Euler). It is also possible to increase stability by going partially implicit (i.e., more implicit) by handling certain quantities implicitly, while treating others explicitly.

Solving the problem implicitly means that you have to solve a set of non-linear equations at each time step. The easiest method of solving a problem like this implicitly is to use the so-called "method of lines" in conjunction with a stiff ODE software package based on the C. W. Gear method. Such a package transparently solves the set on non-linear equations at each time step. The method of lines changes the non-linear PDE into a non-linear set of coupled ODEs, by finite differencing the spatial derivatives, but not the time derivative. So there is an ODE for each grid point. The stiff package automatically integrates these stiff ODEs.

A partially implicit method can be developed for your problem that requires a tri-diagonal set of linear equations to be solved at each time step. This can be done using a tri-diagonal matrix inversion routine. However, there will be some quantities in the difference equations that will have to be handled explicitly in this approach. Still, in my judgement, it would improve stability very substantially (and allow much larger time steps).

The first thing to do is to see whether the existing scheme gives you stability when small (but acceptable) time steps are used.

Chet

 

[member 428835]

Thanks a ton! So I just finished writing the code. It's at the end of the attached pdf. I put plots in the pdf for a time value of 14 (line 8 in the code) and it ran well. However, if I set $tf = 15$ I get some errors, specifically imaginary parts. My guess is from line 37/38, when the square root is introduced. I should say that if I change line 9 to $N = (tf/10)*M^2;$ I can run for $tf$ values of 100 and more.

What is your take on this technique? Is this acceptable or should I go for an implicit scheme?

Also, for your curiosity, the percent change in volume was less than 1.5% for the code and pictures I posted. Changing line 9 as suggested and running for $tf=100$ produces a volume change of about 0.003%.

 

[Chestermiller]

Thanks a ton! So I just finished writing the code. It's at the end of the attached pdf. I put plots in the pdf for a time value of 14 (line 8 in the code) and it ran well. However, if I set $tf = 15$ I get some errors, specifically imaginary parts. My guess is from line 37/38, when the square root is introduced. I should say that if I change line 9 to $N = (tf/10)*M^2;$ I can run for $tf$ values of 100 and more.

What is your take on this technique? Is this acceptable or should I go for an implicit scheme?

Also, for your curiosity, the percent change in volume was less than 1.5% for the code and pictures I posted. Changing line 9 as suggested and running for $tf=100$ produces a volume change of about 0.003%.

I would say it looks pretty OK. I haven't checked the coding, but the solution looks like you expected, correct?

Regarding the need to go implicit, my question is "is computation time a problem?" If not, then I see no need to.

How come you are using such a large M. I would use a value of 20, and then check the accuracy of the solution by redoing the calculation at M = 40. But, I think 20 will be adequate. That should help tremendously with your stability.

What does L(t) look like?

Chet

 

[member 428835]

I haven't checked the coding, but the solution looks like you expected, correct?

Yes, the solution looks good.

Regarding the need to go implicit, my question is "is computation time a problem?" If not, then I see no need to.

Computation time is not an issue.

How come you are using such a large M. I would use a value of 20, and then check the accuracy of the solution by redoing the calculation at M = 40. But, I think 20 will be adequate. That should help tremendously with your stability.

I wanted a large M to reduce error. Having M=20 yields a final $L$ of about 8.4. Changing M=40 yields a final $L$ of about 8.75.

What does L(t) look like?
Chet

I'll attach a pdf.

When changing M=20 and N=M^2 the maximum final run time before crashing is only 20. If I go to M=21 I get complex y values. However, if I let $N=(tf/10)*M$ I can go further in time. What are your thoughts here?

 

[Chestermiller]

I'm a little confused over N. My understanding is that, when you increase N, it makes Δt small. This should increase stability. Yet when you switched from N=M² to $N=(t_f/10)M^2$, you basically reduced N by a factor of about 0.7 (tf went to 100 from 14). This should have made the calculation less stable instead of more stable. Also, when you went to smaller values of M, it should have increased stability at the same value of tf. To make this all easier to follow, you should calculate the value of $\frac{t_fM^2}{(N-1)L^2}$ for each case, and use this to try to map out the region where the solution is stable.

Chet

 

[member 428835]

Let me clarify. For $N=M^2$ and $tf=19$ the code runs well when $M=20$. However, letting $M=40$, all else constant, crashes the code. However, for $N=10*M^2$ and $M=40$ at $tf$ can reach values of 150 without crashing the code. Volume changes at this level are 9.4%.

To make this all easier to follow, you should calculate the value of $\frac{t_fM^2}{(N-1)L^2}$ for each case, and use this to try to map out the region where the solution is stable.
Chet

Why do you choose to use this value?

 

[Chestermiller]

Why do you choose to use this value?

This is Δt/(Δz)² which, aside from the diffusivity factor, is the standard stability parameter for a parabolic equation.

Chet

 

[member 428835]

Thanks. So I checked and the code is roughly stable for $\Delta t / \Delta x^2 < .6$. Additionally, $\Delta t / \Delta x^2$ values of about .65 and up crash. What do you think about this?

As an aside, could you show me where are we actually conserving volume in the pdf? I ask because I was also going to try the same differential equation and the same tip boundary condition but instead of conserving volume apply the following integral constraint $$
-2 \int_0^{tf} \left.h^2 h_z\right|_{z=0} \,d t =\int_0^{L(tf)} \left. h^2 \right|_{t=tf} \, d z$$ I can discretize this for the new $L$. What do you think? Won't this boundary condition affect the previous $L$ derived?

 

[Chestermiller]

The way we have it set up right now, it doesn't automatically conserve volume. But there is a way of setting it up so that it does. Let's start with the differential equation:
$$\frac{\partial h^2}{\partial t}=Z\left(\frac{\partial h^2}{\partial Z}\right)\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{\partial ^2h^3}{\partial Z^2}$$
If we multiply both sides by L, and re-express the first term on the right hand side using integration by parts, we obtain:
$$L\frac{\partial h^2}{\partial t}=\left[\frac{\partial (Zh^2)}{\partial Z}-h^2\right]\frac{dL}{dt}+\frac{2}{3L}\frac{\partial ^2h^3}{\partial Z^2}$$
If we next mover the second term in brackets to the left hand side of the equation, we obtain:
$$\frac{\partial (Lh^2)}{\partial t}=\frac{\partial (Zh^2)}{\partial Z}\frac{dL}{dt}+\frac{2}{3L}\frac{\partial ^2h^3}{\partial Z^2}$$
The finite difference form of this (that will automatically conserve volume) is:

$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}$
This yields:

$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]Δt$

If I did this right, it should give automatic conservation of volume.

Chet

 

[member 428835]

If I did this right, it should give automatic conservation of volume.

I follow your math, but how does it now conserve volume? It seems like you have only written it differently, but have not introduced anything new. Am I missing something?

 

[Chestermiller]

I follow your math, but how does it now conserve volume? It seems like you have only written it differently, but have not introduced anything new. Am I missing something?

In my next-to-last equation, if you add the equations for all the grid points together, the right hand sides will cancel out. This will give you that the trapazoidal rule integration of h^2LdZ for the volume will have a time derivative of zero.

Chet

 

[member 428835]

In my next-to-last equation, if you add the equations for all the grid points together, the right hand sides will cancel out.
Chet

This is not immediately obvious to me. I'm happy to crunch the numbers, but I can't see the cancelation.

This will give you that the trapazoidal rule integration of h^2LdZ for the volume will have a time derivative of zero.
Chet

But h^2LdZ is not volume. Volume was h^2dZ.

 

[Chestermiller]

But h^2LdZ is not volume. Volume was h^2dZ.

No way. h^2dz=h^2LdZ is volume.

Chet

 

[Chestermiller]

This is not immediately obvious to me. I'm happy to crunch the numbers, but I can't see the cancelation.

It works kinda like this: (A-0) + (B-A) + (C-B) + (D-C) + (0 - D) = 0

 

[member 428835]

No way. h^2dz=h^2LdZ is volume.

Chet

Aaaaaaand now I feel stupid. Thanks for clarifying; this makes perfect sense now!

 

[member 428835]

$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}$
This yields:
$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]Δt$

If I did this right, it should give automatic conservation of volume.

I think there is a mistake. I'll begin where you did
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{L Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]\frac{ΔtL}{L(t+Δt)}$

I checked back through the works but can't figure out if the $L$'s are evaluated at $t$ or $t+\Delta t$. I left it as simply $L$ since I am not sure. What do you think?

I now have the reflective boundary condition as $$y_{new}(1) = \frac{L_{old}}{L_{new}}y(1)+\left[ \frac{y(2)}{2 L^2}\frac{dL^2}{dt}-\frac{4 y^{3/2}(1)}{3 L^2 \Delta z^2}\right] \frac{L \Delta t}{L_{new}}$$ where 1 implies the first element of y. Do you agree with this? I should add that everywhere listed as simply $L$ I ran in the code as $L(t+\Delta t)$ yet volume isn't conserved. Any ideas?

 

[Chestermiller]

I think there is a mistake. I'll begin where you did
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{dL}{dt}+\frac{2}{3L}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$\frac{L(t+Δt)h^2(t+Δt,Z)-L(t)h^2(t,Z)}{L Δt}=\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{2ΔZ}\frac{1}{L}\frac{dL}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2} \implies$
$h^2(t+Δt,Z)=\frac{L(t)}{L(t+Δt)}h^2(t,Z)+[\frac{(Z+ΔZ)h^2(t,Z+ΔZ)-(Z-ΔZ)h^2(t,Z-ΔZ)}{4ΔZ}\frac{1}{L^2}\frac{dL^2}{dt}+\frac{2}{3L^2}\frac{h^3(t,Z-ΔZ)-2h^3(t,Z)+h^3(t,Z+ΔZ)}{(ΔZ)^2}]\frac{ΔtL}{L(t+Δt)}$

I checked back through the works but can't figure out if the $L$'s are evaluated at $t$ or $t+\Delta t$. I left it as simply $L$ since I am not sure. What do you think?

I now have the reflective boundary condition as $$y_{new}(1) = \frac{L_{old}}{L_{new}}y(1)+\left[ \frac{y(2)}{2 L^2}\frac{dL^2}{dt}-\frac{4 y^{3/2}(1)}{3 L^2 \Delta z^2}\right] \frac{L \Delta t}{L_{new}}$$ where 1 implies the first element of y. Do you agree with this? I should add that everywhere listed as simply $L$ I ran in the code as $L(t+\Delta t)$ yet volume isn't conserved. Any ideas?

Yes.

It doesn't matter what we do with the L's, since each of the spatial terms on the right hand side have to integrate to zero individually. The reason that they didn't in your calculation is that I used a finite difference approximation in the first term on the right hand side doesn't integrate to zero (my mistake). The approximation I should have used was:
$$\frac{\partial (Zh^2)}{\partial Z}=\frac{(Z+\frac{ΔZ}{2})\frac{(h^2(Z+ΔZ)+h^2(Z))}{2}-(Z-\frac{ΔZ}{2})\frac{(h^2(Z-ΔZ)+h^2(Z))}{2}}{ΔZ}=\frac{(Z+\frac{ΔZ}{2})(h^2(Z+ΔZ)+h^2(Z))-(Z-\frac{ΔZ}{2})(h^2(Z-ΔZ)+h^2(Z))}{2ΔZ}$$

If you try it with this, I guarantee that volume will be conserved.

Chet

 

[member 428835]

The approximation I should have used was:
$$\frac{\partial (Zh^2)}{\partial Z}=\frac{(Z+\frac{ΔZ}{2})\frac{(h^2(Z+ΔZ)+h^2(Z))}{2}-(Z-\frac{ΔZ}{2})\frac{(h^2(Z-ΔZ)+h^2(Z))}{2}}{ΔZ}=\frac{(Z+\frac{ΔZ}{2})(h^2(Z+ΔZ)+h^2(Z))-(Z-\frac{ΔZ}{2})(h^2(Z-ΔZ)+h^2(Z))}{2ΔZ}$$

If you try it with this, I guarantee that volume will be conserved.

Chet

I've redone the code with this expression and no luck. In fact, now the slope is not zero at Z=0. I've attached the code so you can see for yourself. Am I doing something terribly wrong here?

 

[Chestermiller]

I've redone the code with this expression and no luck. In fact, now the slope is not zero at Z=0. I've attached the code so you can see for yourself. Am I doing something terribly wrong here?

I'm going to look over the code and get back with you.

Chet

 

[Chestermiller]

Line 43 looks incorrect to me. I get:

ynew( 1 ) = Lold/Lnew*y ( 1 )+((y ( 2 )+y ( 1 ) ) /(4Lnew^2) * DL2+4/3*(y3( 2 )-y3 ( 1 )) /(Lnew^2 dz ^2) ) dt

My recommended corrections are in bold and italic.
Regarding line 44 (and maybe other lines), I definitely would not overwrite any of the y's until all the ynew's are calculated.

Regarding line 23, this is an incorrect application of the trapazoidal rule. The equation should read:

iv = iv+((h^2( i )+h^2( i+1) ) /2) dz

That's about all I checked so far.

Chet

 

[member 428835]

Line 43 looks incorrect to me. I get:

ynew( 1 ) = Lold/Lnew*y ( 1 )+((y ( 2 )+y ( 1 ) ) /(4Lnew^2) * DL2+4/3*(y3( 2 )-y3 ( 1 )) /(Lnew^2 dz ^2) ) dt

My recommended corrections are in bold and italic.
Regarding line 44 (and maybe other lines), I definitely would not overwrite any of the y's until all the ynew's are calculated.

Regarding line 23, this is an incorrect application of the trapazoidal rule. The equation should read:

iv = iv+((h^2( i )+h^2( i+1) ) /2) dz

That's about all I checked so far.

Chet

Awesome! Not sure how I missed this but it looks great now, and conserves volume! You've been a huge help, Chet! I'll post a final pdf when I finish it so you can see how it all looks!

 

[member 428835]

Hey Chet, I have a quick question. When transforming from $Z$ to $z$ I did the following: $Z=L(tf)z$, but now I'm not sure why I used the last value of $L$. Should $L$ be evaluated at the final time? Could you explain your reasoning?