Finite Differencing Error Help

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

 

[Chestermiller]

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?

You did it so you could always have a grid point exactly at the moving boundary, making it convenient to apply any boundary condition there. Also, you always have the same number of grid points, and they are uniformly distributed between the boundaries. You've taken the complexity out of applying the moving boundary condition and transferred it to the differential equation where it can be handled much more easily.

Chet

 

[member 428835]

Thanks! For some reason I was wondering which L I should use to transform back from Z to z but obviously it's the most recent L. Thanks Chet! You're definitely brilliant!

 

[member 428835]

Hey Chet, sorry to bother you again about this problem, but how much more effort would need to be put into the finite difference scheme to go implicit? I've never done this before and am curious about it.

 

[Chestermiller]

Hey Chet, sorry to bother you again about this problem, but how much more effort would need to be put into the finite difference scheme to go implicit? I've never done this before and am curious about it.

It depends. What you would do would be to only discretize the differential equations in the spatial direction, but not in the time direction. So then you would have a coupled set of non-linear first order ODEs in time for the y values at the M grid points (plus an additional first order ODE for L²). This approach is called the Method of Lines. To solve these equations implicitly, you would typically need a stiff ODE equation solver subroutine. The problem would be done in FORTRAN, and you would be using a commercially available stiff package like the ones available on IMSL or online (e.g., DASSL). Or, if Matlab has a stiff equation solver, you could use that.

Chet

 

[member 428835]

Thanks a ton Chet! I really appreciate your help.

 

[Strum]

Is #86 the newest code? I would like to look it over when it is fully done :)

 

[member 428835]

after chet's correction, yep!